Application of GRAAL model to the resumption of International Simple Glass alteration

Supplementary Figure 6 :Supplementary Figure 7 :

 67 Application of GRAAL model to the resumption of International Simple Glass alteration M . F o u r n i e r , P . F r u g i e r , S . G i n P age 3/7 Supplementary Figure 5: Comparisons between modeled and experimental concentrations at pH 10.4. Experimental B, Si, Na, and Al 1 (diamonds) and modeled (dashed line for  = 1 and dotted line for   0) concentrations in (a) unseeded and (b) seeded tests conducted at pH 10.4 maintained by adding NaOH in static conditions, at 90°C, S/V = 1,770 m -1 . Application of GRAAL model to the resumption of International Simple Glass alteration M . F o u r n i e r , P . F r u g i e r , S . G i n P age 4/7 Comparisons between modeled and experimental concentrations at pH 11. Experimental B, Si, Na, and Al 1 (diamonds) and modeled (dashed line for  = 1 and dotted line for   0) concentrations in (a) unseeded and (b) seeded tests conducted at pH 11 maintained by adding NaOH in static conditions, at 90°C, S/V = 1,770 m -1 . Application of GRAAL model to the resumption of International Simple Glass alteration M . F o u r n i e r , P . F r u g i e r , S . G i n P age 5/7 Comparisons between modeled and experimental concentrations at pH 11.3. Experimental B, Si, Na, and Al 1 (diamonds) and modeled (dashed line for  = 1 and dotted line for   0) concentrations in (a) unseeded and (b) seeded tests conducted at pH 11.3 maintained by adding NaOH in static conditions, at 90°C, S/V = 1,770 m -1 .

SUPPLEMENTARY INFORMATION -Application of GRAAL model to the resumption of International Simple Glass alteration M. Fournier, P. Frugier, S. Gin P age 1/7

Application of GRAAL model to the resumption of International Simple Glass alteration

Maxime Fournier, Pierre Frugier, Stéphane Gin

SUPPLEMENTARY INFORMATION

Supplementary Note 1 Assumptions for modeling zeolites nucleation-growth Supplementary Figure 1: In the proposed approach, nucleation occurs at a constant rate on the gel surface area (glass surface area). The nuclei surface area available for zeolites growth is therefore proportional to the time t:

The zeolites growth occurs is in a single direction, forming needles.

Supplementary Figure 2:

In seeded tests, the growth surface area Sc, corresponding to the seeds surface area introduced into the medium at the initial time, is independent of time. The seeds growth is represented by the unidirectional growth of a set of needles.

Supplementary Note 2 Gel description in an alkaline environment

Supplementary Figure 3: In order to integrate the RA phenomenon to the GRAAL model, a new gel composition domain was defined. This domain integrates a passivating (PRI) and three non-passivating endmembers (SiAlCa, SiZrNa, and SiZr0.1Na0.2). Two secondary phases (Na-P2 zeolites and CSH0.8) also play a role in controlling the activities in solution.

Phase controlling activity during

Supplementary Discussion

The solubility of the Na-P2 zeolite seeds was evaluated by placing synthesized seeds 1 in NaOH solutions of different molarities for one month at 90 °C. The ratio between the seeds "geometric" surface area (0.25 m 2 •g -1 ) and the solution volume was approximately 2•10 -3 m -1 . The concentrations measured in solution were fairly constant, indicating the early onset of saturation conditions in solution (Supplementary Figure 8). Even at the lowest pH tested (therefore the farthest pH from the synthesis conditions), the dissolution of the seeds was low and the stationary state was reached quickly. The altered zeolite thickness increased with the increase in the pH. Moreover, the Si/Al ratio in solution increase, thus deviating from 1.7, the ratio of the seeds (Supplementary Table 2). The seeds dissolution was not congruent, implying the slow formation of a phase with a stoichiometry different from the starting mineral. Since the maximum equivalent thickness of Si was only 300 nm after one month of leaching, the quantities formed of this (or these) phase(s) remained low, explaining that their identification and characterization could not be completed (not visible by SEM after 330 days at pH ≈ 9). The solubility constant of the Na-P2 zeolite was calculated using the CHESS code for different pH values. This calculation required the implementation of the following dissolutionprecipitation reaction in the database, using a stoichiometry measured elsewhere 1 :

For all the pH values studied, we calculated a solubility constant value at each sampling time (1, 3, 7, 14, and 28 days respectively). The results obtained were similar for the same pH but varied substantially from one pH to another (Supplementary Figure 9). reasonably describes the equilibrium between the Na-P2 zeolite and the NaOH solutions (Supplementary Figure 10). This log K value was used in the calculations presented in the main text.


Keywords: Neutronics, Transport, 2D/1D method, Homogenization, Energy condensation, Sodium Fast Reactor

Fast resolution of the Boltzmann transport equation over a nuclear reactor core presupposes the definition of homogenized and energy collapsed cross sections. In modern sodium fast reactors that rely on heterogeneous core designs, anisotropy in the neutrons propagation cannot be neglected so three-dimensional models should be preferred to compute those effective cross sections. In this paper, the 2D/1D approximation is used to avoid computationally expensive 3D calculations while preserving consistent angular representations of the neutron flux. An iterative procedure is defined to solve the 2D/1D equations and produce coarse group homogenized cross sections that account for 3D transport effects. Accuracy of the algorithm is tested on a realistic model of the ASTRID core showing very good results against Monte Carlo simulations for all neutronic parameters (eigenvalue, sodium void worth and fission map distribution).

I. INTRODUCTION

Next generation of nuclear reactors is asked to meet very high safety criteria in order to prevent as much risk of accident as possible. In the field of neutronics, reaching those standards presupposes the demonstration of a good understanding of the physical processes that take place in the reactor core. Given the progress that has been made in computer science, numerical simulation nowadays provides an efficient tool to perform such a demonstration and is therefore widely used.

Consequently, the research community has concentrated many efforts in providing detailed and precise solutions to three-dimensional whole-core transport problems. While taking a significant step forward in the quest for accuracy, those efforts also showed that obtaining such solutions often required long hours of calculation over large-scale parallel systems with hundreds and even thousands of processors. As a result, they are not well adapted to the current industry's needs which demand a large number of calculations in a limited amount of time and with limited computer resources [START_REF] Sanchez | Prospects in deterministic three-dimensional whole-core transport calculations[END_REF]. For practical application then, we may assess with confidence that full core calculations will still be performed for some time over coarse spatial and energy meshes. They will therefore rely on homogenized multigroup effective cross sections.

Those effective parameters are classically obtained from two-dimensional lattice calculations in which neutron fluxes are computed over axially infinite and radially periodic assembly patterns and used to weight the input cross sections. Taking advantage of nearly periodic core designs, limited axial heterogeneities and neutron' small mean free paths, the lattice paradigm has proven its efficiency for pressurized water reactors (PWR) applications.

On the contrary, current tendency for modern sodium-cooled fast reactors (SFR) is to present strong axial heterogeneities. While allowing a significant reduction of the reactivity void worth compared to homogeneous designs, this feature however questions the validity of the above mentioned procedure for cross sections generation. In particular, the quite large mean free paths of neutrons in sodium-cooled technologies gives birth to particle exchanges between distinct axial layers and therefore produces spectrum shifts along the vertical axis and anisotropic structures in the neutron flux.

More precisely, it has been shown in a previous work [START_REF] Faure | Neutronic calculation of an axially heterogeneous ASTRID fuel assembly with APOLLO3 R : Analysis of biases and foreseen improvements[END_REF] that large errors could indeed be found when two-dimensional lattice parameters were used to compute neutron fluxes over axially heterogeneous SFR cores such as the French CFV (coeur a faible vidange or low sodium void effect core) prototype ASTRID. The conclusion of the analysis was that infinite lattice calculations were unable to yield consistent effective cross sections because they neglect axial modes of the angular flux. Moreover, it was proven that biases could be significantly reduced if such modes were incorporated in the calculation and used in the cross sections weighting process.

A straightforward way to access full angular modes of the flux that are representative of the core situation is the application of the method of characteristics (MOC) to three-dimensional transport problems [START_REF] Santandrea | Accelerated polynomial axial expansions for full 3D neutron transport MOC in the APOLLO3 code system as applied to the ASTRID fast breeder reactor[END_REF]. However, considerable efforts are needed to solve the 3D MOC equations, making its use rather inappropriate in the scope of cross sections generation were fluxes are used only as weight functions.

In this paper, consistent information for cross sections generation is retrieved from a 2D/1D approximation to the 3D transport equation. Based upon the original 2D/1D fusion method [START_REF] Cho | A fusion technique of 2D/1D methods for threedimensional whole-core transport calculations[END_REF], its basic idea is to couple several two-dimensional lattice calculations through axial leakage computed over one-dimensional models. Taking advantage of the angular information that is contained in the axial leakage, it is shown that three-dimensional transport effects can be successfully stored in effective cross sections.

The aim of the paper is to present this 2D/1D cross sections weighting algorithm and to stress its accuracy in realistic situations. Its organization is as follows: Sec. II addresses the theory of effective cross sections generation while Sec. III presents the 2D/1D algorithm together with its practical implementation in the APOLLO3 R code [START_REF] Schneider | APOLLO3: CEA/DEN deterministic multi-purpose code for reactor physics analysis[END_REF]. In Sec. IV, numerical results are presented and discussed. The paper ends with a general conclusion and perspectives for future work.

II. EFFECTIVE CROSS SECTIONS GENERATION

In this section, the transport equation is introduced and the theory of cross sections averaging is discussed.

II.A. The transport equation

The steady state Boltzmann transport equation for the neutron flux ψ reads

Ω • ∇ + Σ ψ = q (1)
where Σ is the total macroscopic cross section and q a source term. Equation (1) stands in a geometrical domain D, for all neutron directions and at all energies (i.e. ( r, E, Ω) ∈ D × R + × S 4π ) provided that the incoming neutron density at the boundary is known (i.e. ψ = ψ in for ( r, Ω) ∈ ∂D × S - 2π ).

In a reactor core, fission and scattering reactions are the main contributors to the neutron source. Introducing the scattering H and production F operators q reads

q = H + 1 k F ψ (2)
In that case, equation ( 1) is casted in the form of an eigenvalue problem. Assuming k is the largest positive eigenvalue, the eigen-pair (k, ψ) is called fundamental mode. It is the solution that is handled by most numerical methods and the one we are interested in.

For practical applications, the energy interval is split into N g bins R + = g E g ; E g+1 and continuous energy data are replaced by their group averaged value. Flux ψ is sought under the form of a vector ψ = (ψ g ) g∈ 1;Ng .

In this well known multigroup formalism, the scattering operator H is written in terms of the Legendre moments of the scattering cross section (Σ g →g s,l

) l∈N as

Hψ g ( r, Ω) = Ng g =1 +∞ l=0 2l + 1 4π Σ g →g s,l ( r)
+l m=-l

φ g lm ( r)R lm ( Ω) (3)
The production operator F is also developed in terms of a production cross section νΣ g →g f as

F ψ g ( r, Ω) = 1 4π Ng g =1 νΣ g →g f ( r)φ g 00 ( r) (4)
Here νΣ g →g f is a non-standard but compact notation that will be used throughout this paper for the sake of simplicity. It stands for the sum over fissile isotopes of the product of fission spectra χ times multiplicity ν times fission cross section Σ f .

Real spherical harmonics (R lm ) l∈N m∈ -l;l have also been used as a projection basis for angular dependent functions. Angular flux moments are defined as

φ g lm ( r) = S4π d 2 Ω ψ g ( r, Ω)R lm ( Ω) (5)

II.B. Cross sections averaging

A typical SFR fuel assembly is composed of 250 fuel pins requiring thus at least 600 radial volumes to discretize fuel, clad, coolant and internal structures. Considering that there are approximately 500 assemblies in the core and that 100 distinct axial layers are required to catch the axial variations of the flux, a total of 30 million spatial volumes is a minimum. If we also suppose that 250 angles are required for the Ω variable and that N g = 2000 is a realistic number of energy groups (when resonances are accounted for with a self-shielding formalism), we get the digit of 1.5 10 13 unknowns for the flux. Assuming a simple precision floating-point storage, the latter represents a total of 60 terabytes of memory to which must be appended the weight of additional data such as cross sections, geometry description or transport acceleration factors.

Given that such a huge amount of computer resources is not suited for industrial applications, cross sections averaging is required to reduce the number of energy groups (energy condensation)

and spatial mesh cells (homogenization) in whole core calculations. Its basic idea is to preserve reaction and leakage rates between detailed and coarse models of the same transport problem.

In order to derive analytical expressions for effective cross sections, we first introduce compact notations for spatial integration and group summation

< f g > iG = Di d 3 r g∈G f g ( r, Ω) (6)
< f g →g > iG G = Di d 3 r g ∈G g∈G f g →g ( r, Ω) (7)
where D i ⊂ D is a generic sub-domain and G = {g 1 , . . . , g p }, G = {g 1 , . . . , g p } are energy group indexes sets. We then define the following quantities

Σ G i ( Ω) = < Σ g ψ g > iG < ψ g > iG (8) Σ G →G s,lm,i = < Σ g →g s,l φ g lm > iG G < φ g lm > iG (9) νΣ G →G f,i = < νΣ g →g f φ g 00 > iG G < φ g 00 > iG (10)
At first sight, such definitions are good candidates for effective total, scattering and production cross sections. They however suffer from two well known limitations that are here recalled:

1. They generally do not preserve exactly the neutronics balance because < ψ g > iG differs from the integral of the coarse-mesh transport solution. Two exceptions might be found, when only energy collapsing is considered [START_REF] Vidal | An improved energy-collapsing method for core-reflector modelization in SFR core calculations using the PARIS platform[END_REF] or when a whole space homogenization is performed with zero net leakage boundary conditions [START_REF] Sanchez | Assembly homogenization techniques for core calculations[END_REF]. In all other cases, preservation of leakage rates is not guaranteed.

2. Averaging equation ( 8) with the angle dependent function ψ introduces anisotropy in the total cross section. As a result, Σ G i does not define a transport equation over an isotropic medium and is thus not compatible with most neutron transport codes. Similarly, the dependence of the scattering cross section with azimuthal index m in equation ( 9) is not standard.

To address the first issue, SPH factors have been proposed in the past by Hebert [START_REF] Hébert | Development of a third-generation superhomogeneisation method for the homogenization of a pressurized water reactor assembly[END_REF] for partial homogenization and pin-by-pin power reconstruction. In the same time, Smith also introduced current (and latter flux) discontinuity factors [START_REF] Smith | Assembly homogenization techniques for light water reactor analysis[END_REF] to deal with general boundary conditions (black box homogenization).

As for the second point, a standard solution is the flux-volume collapsing technique that consists in replacing angular dependent weights ψ and φ lm by the isotropic scalar flux φ 00 . If this approximation clearly implies the loss of any angular information, consistency can be recovered in the framework of diffusion theory because diffusion coefficients can be constructed in such a way as to preserve linearly anisotropic transport effects. When transport theory is used for whole core calculations however, high order angular effects should be stored in effective cross sections.

In the generalized energy condensation theory introduced by Rahnema and co-workers [START_REF] Rahnema | Generalized energy condensation theory[END_REF][START_REF] Douglass | Consistent generalized energy condensation theory[END_REF], the angular dependence of the total cross section is shifted into a variation term and transferred to the right hand side of the transport equation. Following the same idea, Vidal incorporated this variation term into the within group scattering rate [START_REF] Vidal | An improved energy-collapsing method for core-reflector modelization in SFR core calculations using the PARIS platform[END_REF][START_REF] Vidal | APOLLO3 Homogenization Techniques for Transport Core Calculations -Application to the ASTRID CFV Core[END_REF] and used a least-square minimization technique to define a conventional scattering cross section.

In this paper, we will follow this latter approach and take the following definitions for effective total and l-th order scattering cross sections:

Σ G i = < Σ g φ g 00 > iG < φ g 00 > iG (11) Σ G →G s,l,i = +l m=-l Σ G →G s,lm,i < φ g lm > iG 2 +l m=-l < φ g lm > iG 2 + δ G,G Σ G i -Σ G l,i (12)
where δ G,G is the Kronecker delta, Σ G →G s,lm,i is taken from equation ( 9) and Σ G l,i reads

Σ G l,i = +l m=-l < φ g lm > iG < Σ g φ g lm > iG +l m=-l (< φ g lm > iG ) 2 (13)
Equation ( 10) defines a standard cross section and is therefore used for the production probability.

In the following, the combination of equations ( 10), ( 11) and ( 12) will be referred to as flux moments weighting. It allows a better conservation of anisotropic transport effects and is therefore expected to improve the standard (scalar) flux-volume weighting technique.

III. THE 2D/1D METHOD FOR CROSS SECTIONS GENERATION

The 2D/1D terminology refers to a class of computational methods that proposes to solve the three-dimensional Boltzmann transport equation by resorting to transverse integration in order to reduce the phase-space dimensions. In this section, we present how the 2D/1D method can be used in the scope of effective cross sections generation.

The original idea comes from Cho and colleagues that developed a fusion method for the CRX code coupling two-dimensional radial MOC calculations one to another through axial leakage retrieved from one-dimensional S N solutions [START_REF] Cho | A fusion technique of 2D/1D methods for threedimensional whole-core transport calculations[END_REF]. At the same time, a similar method was also implemented in the DeCART code by another group of researchers that proposed to embed the axial leakage calculation into three-dimensional coarse mesh finite difference simulations [START_REF] Cho | Axial SPN and radial MOC coupled whole core transport calculation[END_REF].

In the past few years, a growing interest was found all around the world for the 2D/1D method [START_REF] Ban | Code-to-code comparisons on spatial solution capabilities and performances between nTRACER and the standalone IDT solver of APOLLO3[END_REF][START_REF] Tang | Efficient procedure for radial MOC and axial SN coupled 3D neutron transport calculation[END_REF]. In particular, a mathematical analysis of the 2D/1D equations and their numerical counter part was provided by the MPACT development team [START_REF] Kelley | A consistent 2D/1D approximation to the 3D neutron transport equation[END_REF] and some of the main limitations, such as the negative source issue, were recently lifted [START_REF] Collins | Stability and accuracy of 3D neutron transport simulations using the 2D/1D method in MPACT[END_REF].

A particular feature of the 2D/1D method is that it requires a homogenization procedure to produce the one-dimensional calculation mesh and a de-homogenization method to build radial shapes for axial leakage. In order to get rid of those approximations, some authors proposed to equalize the number of axial calculations to the number of two-dimensional radial cells [START_REF] Lee | 2D/1D fusion method solutions of the three-dimensional transport OECD benchmark problem C5G7 MOX[END_REF].

However, it is our opinion that computational advantages can be taken from the homogenization step so this latter approach will not be considered in this paper.

In the following, the basic theory of the 2D/1D method is presented and choices made for their practical implementation in the APOLLO3 R code are discussed.

III.A. General equations

The starting point is the 3D transport equation ( 1) whose solution is written (k, ψ). For the sake of clarity, energy group indexes are dropped out in this section.

The geometrical domain D is partitioned into axial Z i and radial D r sub-domains such as

D = i Z i × r D r .
The axial partition is chosen so that cross sections are axially invariant in a given layer

Z i i.e. ∀z ∈ Z i , Σ( r) = Σ(x, y) (14)
Equation ( 1) is then successively integrated axially over Z i and radially over D r yielding the following set of equations

Ω • ∇ xy + Σ ψ i = q i -L i ( 15
) (µ∂ z + Σ r ) ψ r = q r -L r (16)
where ψ i , ψ r are integrated fluxes

ψ i (x, y, Ω) = Zi dz ψ( r, Ω) (17)
ψ r (z, Ω) = Dr dxdy ψ( r, Ω) (18)
and L i , L r transverse leakage

L i (x, y, Ω) = µ Zi dz ∂ z ψ( r, Ω) (19) L r (z, Ω) = Dr dxdy Ω • ∇ xy ψ( r, Ω) (20)
∂ z symbolizes the partial derivative towards z and µ is the polar cosine.

Applying operators H and F of equation ( 2) to flux ψ i (resp. ψ r ), one finds the analytical expression for source q i (resp. q r ). Cross section Σ is the one from the original 3D problem because of assumption [START_REF] Ban | Code-to-code comparisons on spatial solution capabilities and performances between nTRACER and the standalone IDT solver of APOLLO3[END_REF]. On the contrary, the definition of Σ r requires a 3D flux homogenization

Σ r = < Σψ > Dr < ψ > Dr (21)
Similar expressions also stand for cross sections that are implicit in q r (scattering and production).

At this point, the above equations suffer from no approximation except assumption [START_REF] Ban | Code-to-code comparisons on spatial solution capabilities and performances between nTRACER and the standalone IDT solver of APOLLO3[END_REF].

They however cannot be avoided if the dependence on the 3D flux ψ is to be dropped out in transverse leakage L i , L r and in the averaged cross section Σ r .

III.B. Approximations for SFRs

This paper addresses the topic of cross sections averaging so geometrical domain D only needs to be a subdomain of the core suited for computation of representative weight functions (flux moments). In the following we therefore suppose that D is a radially reflected 3D pattern (a fuel assembly for instance).

In SFRs neutrons have quite large mean free paths (up to a few centimeters) and it is possible to assume that neutron fluxes are locally flat. We therefore consider no radial partitioning of domain D but only an axial splitting i.e.

D = i Z i × D r .
Under those circumstances, the zero net radial leakage boundary condition (radial reflection) over ∂D r leads to

L r = 0 ( 22
)
If we also neglect the radial dependence of leakage L i and integrate equation ( 19) over D r , we get

L i = µ A xy Zi dz ∂ z ψ r = µ A xy ψ r | zi+ zi- (23)
with A xy = Dr dxdy and

Z i = [z i -; z i +].
Finally, assuming that the change in the radial flux shape is small along Z i , cross section Σ r can be averaged with 2D flux ψ i instead of 3D flux ψ i.e.

Σ r,i = < Σψ i > Dr < ψ i > Dr (24)
Most 2D/1D methods neglect the angular dependences in the homogenized cross sections and use a scalar flux weighting in the above equation. However, the preservation of anisotropic effects in transport calculations requires higher order methods as it was pointed out in Sec. II.B. In this paper, the flux moments weighting technique (i.e. equations ( 10), ( 11) and ( 12)) is used. Other solutions such as polar angle dependent cross sections [START_REF] Jarret | Characterization of 2D/1D accuracy in MPACT[END_REF] might also be found in the literature but their number is limited.

Introducing those approximations in the 2D/1D system, we find that equation ( 15) becomes a 2D equation for flux ψ i with imposed external source L i while equation ( 16) becomes a standard 1D eigenvalue problem for flux ψ r with homogenized cross section Σ r . They form the closed system of coupled 2D/1D equations.

III.C. Interest for cross sections collapsing

The 2D/1D system of equations defines a natural algorithm for producing effective cross sections that preserve 3D transport effects at a relatively low computational cost. A schematic view of this algorithm is presented in Figure 1. It reads:

1. Initialize leakage L i 2. Compute 2D fluxes ψ i with eq. ( 15)

3. Homogenize cross sections Σ r,i with eq. ( 24)

4. Compute 1D flux with eq. ( 16) and ( 22

L i = L 0 i 2D calculations ∀i, Ω • ∇ xy + Σ ψ i = q i -L i Homogenization ∀i, Σ r,i = <Σψi> Dr <ψi> Dr 1D calculation (µ∂ z + Σ r ) ψ r = q r Converged? No Axial leakage ∀i, L i = µ Axy ψ r | zi+ zi-

Yes

Effective cross sections eqs. ( 10), ( 11), ( 12) Fig. 1. 2D/1D algorithm for effective cross sections generation Convergence of the algorithm can be checked on 2D or 1D fluxes, effective cross sections or, more simply, on the 1D eigenvalue. In Sec. IV, this latter test will be used.

The main interest of the 2D/1D equations for cross sections collapsing is the representation of the angular variable. In standard 2D patterns (lattice paradigm), high order angular modes of the flux are often close to zero because anti-symmetric modes with respect to the (xOy) plane do not exist and reflective boundary conditions also alter radial moments. As a consequence, < φ lm > i is close to 0 for l > 0 and equation ( 12) cannot be used to collapse high order cross sections. On the contrary, the 2D/1D framework involve full angular fluxes. So consistent flux moments can be used to collapse cross sections (if such moments exist in the 3D geometry).

Another interesting feature of the 2D/1D algorithm is that it naturally takes into account spectrum shifts along the (Oz) axis while it is not the case of standard 2D lattice calculations.

Moreover the algorithm provides a way to compute effective cross sections in media that have a low (or null) intrinsic neutron source (e.g. fertile blankets or shields) and avoids thus the construction of ad-hoc sources or cluster geometries.

Finally, the 2D/1D approximation avoids full 3D simulations and offers flexibility in two dimensional simulations that can be performed independently. In addition, axially non-conformal geometries can be treated without any restriction.

III.D. Practical implementation

The previous algorithm has been implemented in the APOLLO3 R code to test its validity.

In their most general form, equations ( 15) and ( 16) involve 3D angular representations of fluxes ψ i and ψ r . However, the trajectory sweeping in standard 2D and 1D solvers is done accordingly to the symmetries present in the geometry. As a consequence some angular moments are automatically set to 0 and the 3D angular information is lost. So particular solvers are required.

III.D.1. 2D solver

A special version of the TDT MOC solver [START_REF] Santandrea | An Integral Multidomain DPN Operator as Acceleration Tool for the Method of Characteristics in Unstructured Meshes[END_REF] was recently developed in APOLLO3 R to solve the heterogeneous B equations in 2D geometries. Those equations come from the theory of neutron leakage in infinite lattices [START_REF] Deniz | The Theory of Neutron Leakage in Reactor Lattices[END_REF] and read

Ω • ∇ + Σ + i B • Ω ψ = q (25)
The TDT-B solver sweeps polar angles in the full [0; π] range so 3D angular representations are allowed for the flux. If the buckling vector B is set to 0, equation ( 25) is formally identical to equation ( 15) so the TDT-B solver can be used to solve the 2D equations.

III.D.2. 1D solver

When approximations of Sec. III.B are made, radial leakage are set to zero in equation [START_REF] Kelley | A consistent 2D/1D approximation to the 3D neutron transport equation[END_REF].

So the right hand side of the equation is not impacted by 3D angular effects. In addition, the moments weighting technique does not introduce particular angular dependences in the homogenized 1D cross sections. As a result, a standard 1D solver can be used. In this paper, the short characteristics IDT solver [START_REF] Zmijarevic | IDT solution to the 3D transport benchmark over a range in parameter space[END_REF] of APOLLO3 R is chosen.

III.D.3. 2D/1D coupling

Axial leakage L i are computed from 1D interface angular fluxes with equation [START_REF] Coste-Delclaux | GALILEE: A nuclear data processing system for transport, depletion and shielding codes[END_REF]. They are then projected over Legendre polynomials P l and identified to rotationally (Oz)-invariant angular moments for 2D calculations

L i,lm =        1 A xy +1 µ=-1 dµ P l (µ)µ [ψ r (z, µ)] zi+ zi-if m = 0 0 otherwise (26)

III.D.4. Fix-up for negative sources

There is no guarantee that the total source in the right hand side of equation ( 15) is positive.

As a result, instabilities and even negative fluxes can be found when axial leakage are strong enough. In reference [START_REF] Collins | Stability and accuracy of 3D neutron transport simulations using the 2D/1D method in MPACT[END_REF], this issue is addressed with transverse leakage splitting and the negative source is turned into an additional absorption term.

In this work however, it has been chosen to avoid intrusions in flux solvers so the positiveness of q i -L i is never tested. Stability of the algorithm is nevertheless guaranteed by neglecting physical leakage i.e.

L i = 0 when -L i < 0 (27)
In sub-critical and non multiplicative 2D layers Z i × D r (shields or fertile blankets), approximation (27) does not affect the stability of the algorithm because the existence of a neutron source q i -L i is ensured for all energy groups: -L i > 0 at high energies (arrival of high energy neutrons from main fissile zones) and q i > 0 at low energies (down-scattering effect).

For over-critical 2D layers (fuel) on the contrary, approximation (27) creates the condition of an over-critical pattern with positive external source. In such a situation, no solution of the transport equation exists so any iterative resolution of equation ( 15) is condemned to diverge. The physical interpretation is that leakage cannot be neglected in over-critical layers. If positiveness of the source q i -L i does not want to be tested, an artifact must then be used to replace real leakage by fictitious ones. The arbitrary choice we made in this paper is to use the heterogeneous B equations ( 25) to account for leakage out of 2D overcritical layers. In other words, real leakage are replaced by a buckling vector B in such patterns. Because B accounts for axial leakage, we choose its orientation along (Oz). Its amplitude B 2 is also adjusted to ensure k(B 2 ) = 1. If those choices are purely arbitrary, they considerably simplify the 2D/1D algorithm since they imply that overcritical layers can be calculated once and for all at the beginning of the iterative strategy.

Reflections are currently under consideration to mitigate those rather crude approximations.

Despite all, we shall see in the next section that quite good results could be obtained.

IV. NUMERICAL RESULTS FOR A SFR ASSEMBLY CALCULATION

IV.A. Problem description

The benchmark chosen for analysis of the 2D/1D algorithm is a hexagonal SFR fuel assembly representative of the ASTRID core. The geometry shows strong axial heterogeneities (CFV design) so axial flux gradients are expected.

The axial layout of the core is depicted in Figure 2.a. C1 stands for fissile (U,Pu)O 2 material while FCA is for fertile UO 2 media. PLN and PNS are French acronyms for sodium plenum and axial neutronic protection respectively. Radial mesh of the fuel pins lattice is presented in Figure 2.b while plenum is shown in Figure 2.c. Dimensions and compositions of the benchmark are available in reference [START_REF] Faure | Neutronic calculation of an axially heterogeneous ASTRID fuel assembly with APOLLO3 R : Analysis of biases and foreseen improvements[END_REF] except for the axial protection (Figure 2.d) that was not modeled is that former work. The latter is composed of a lower plug and a 57.5 cm height column of B 4 C absorber material (7.5 cm 90% 10 B enriched + 50 cm natural boron).

The axial partition of domain D for the 2D/1D algorithm follows the axial layout of the assembly leading to thirteen 2D layers (TDT calculations). In IDT 1D calculations, the axial mesh is refined to ensure a converged flux solution (sub-plane scheme). Taking advantage of the algorithm flexibility, all radial heterogeneities are explicitly described in 2D layers.

Nuclear data are taken from the JEFF-3.1.1 evaluation and prepared into a 1968 group energy structure using the dedicated GALILEE processing tool [START_REF] Coste-Delclaux | GALILEE: A nuclear data processing system for transport, depletion and shielding codes[END_REF]. A P 3 angular order is chosen for the scattering kernel.

Before any flux calculation, cross sections are self-shielded over the exact 2D geometries using the Tone method of APOLLO3 R , which has proven to give very good results for SFR calculations and even comparable with reference subgroup methods [START_REF] Faure | Neutronic calculation of an axially heterogeneous ASTRID fuel assembly with APOLLO3 R : Analysis of biases and foreseen improvements[END_REF][START_REF] Mao | A new Tone's method in APOLLO3 R and its application to fast and thermal reactor calculations[END_REF].

IV.B. Algorithm validation

Before considering the generation of effective cross sections, the behavior of the 2D/1D algorithm itself is tested. In particular, convergence is checked on the 1D eigenvalue k and the influence of initialization is discussed.

IV.B.1. Convergence

Table I shows the reactivity error towards reference value

∆ρ = 10 5 1 k ref - 1 k (28)
as a function of the number of iterations in the 2D/1D algorithm. An iteration, identified in the following with superscript n, is defined as a complete loop in Figure 1 starting from 2D calculations.

For overcritical layers we recall that the B equations ( 25) are solved once and for all in the first iteration. The B leakage rate is used as an initialization source L (0) i for all other planes (FCA, PLN. . . ). 1D cross sections Σ

(n)
r,i are homogenized with 2D flux moments i.e. with equations ( 10), [START_REF] Douglass | Consistent generalized energy condensation theory[END_REF] and [START_REF] Vidal | APOLLO3 Homogenization Techniques for Transport Core Calculations -Application to the ASTRID CFV Core[END_REF].

Reference

k = 1.08183 ± 1 pcm 2D/1D Iteration n o 1 2 3 ∆ρ (pcm) +164 +69 +69 TABLE I Convergence of the 1D eigenvalue in 2D/1D algorithm
We see that the algorithm converges quickly and that the reactivity discrepancy towards reference Monte Carlo results is limited (∆ρ = 69 pcm). More precisely, we observe that only two iterations are necessary before matching convergence: ∆ρ (3) = ∆ρ (2) . The interpretation is that 1D homogenized cross sections Σ r,i are in fact not affected by small variations in the sources L i that are used for computation of 2D homogenization fluxes ψ i . In other words, the first 1D calculation contains enough information to compute accurate-for-homogenization sources L i . Furthermore, it has been observed that non-multiplicative and sub-critical layers flux calculations are not even necessary in first iteration. In fact, cross sections Σ

IV.B.2. Influence of buckling normalization for over-critical layers

As mentioned in Sec. III.D, the 2D/1D equations are not used for overcritical layers but equation ( 25) is solved at the beginning of the iterative process. The buckling B is directed along the (Oz) axis and its modulus | B| adjusted so that k(B 2 ) = 1.

However, this choice is purely arbitrary and it could as well be decided to adjust more subtly | B|, for instance scaling k(B 2 ) to the current iteration 1D eigenvalue . This dynamic scaling has been tested but it was found that its impact on the final result was negligible (less than 5 pcm) compared to the increased number of 2D flux calculations it required.

IV.B.3. Optimization

With the above considerations, a significant number of flux calculations can be avoided without impacting accuracy of the final results. In fact, we are able to define an algorithm in which only one flux calculation is performed for each 2D layer. It reads:

1. Solve B equations (25) for overcritical layers imposing k(B 2 ) = 1 2. Initialize ψ i = 1 for all other 2D layers 3. Homogenize cross sections Σ r,i with eq. ( 24) 4. Compute 1D flux with eq. ( 16) and ( 22

IV.C. Assembly calculation with effective cross sections

In this section, the optimized algorithm is used to produce effective cross sections. The 1968 group energy structure is condensed into 33 groups and homogenization is performed over each axial layer Z i . A library with thirteen sets of 33 group homogenized cross sections is therefore built.

With this library, we are capable of defining a simplified model of the initial benchmark over which the transport equation ( 1) can be solved. This model is chosen to stress the 2D/1D averaged cross sections against reference continuous energy heterogeneous 3D Monte Carlo simulations.

The S N MINARET solver [START_REF] Moller | MINARET, a deterministic neutron transport solver for nuclear core calculations[END_REF] of APOLLO3 R is used to solve the transport equation over the simplified fuel assembly model. Fine spatial and angular discretizations (2.5 cm axial mesh and 144 angle directions) ensure that the flux is converged.

Results for eigenvalues and fission rate distributions, in nominal and voided conditions, are presented in Table III andFigure 3. Equation ( 28) is used for the reactivity discrepancy between APOLLO3 R and TRIPOLI4 R .

To show the need for high order angular moments, results with (scalar) flux volume weighted cross sections are also reported in the second line of Table III. This is done using the same 2D/1D flux distribution but replacing φ lm by φ 00 in equation [START_REF] Vidal | APOLLO3 Homogenization Techniques for Transport Core Calculations -Application to the ASTRID CFV Core[END_REF] It is seen that very good performances are obtained with 2D/1D cross sections when high order angular modes of the flux (moments φ lm ) are used. The results reproduce the Monte Carlo ones within a 60 pcm range for reactivity and within a 1% relative discrepancy for the fission rate distribution. A slightly larger error is found in coolant voided fertile material but in such regions the fission rate is relatively low.

As for the 127 pcm reactivity difference that is found between the 33 group MINARET calculation (∆ρ = -58 pcm) and the 1968 group IDT results (∆ρ = +69 pcm in Sec. IV.B), it mainly accounts for the reduction in the number of energy groups and, to a lesser extent, for differences in flux solvers.

In the other hand, we observe that very large errors arise when only the isotropic component of the flux (φ 00 ) is used to weight cross sections. In particular, a +625 pcm reactivity discrepancy is found in voided conditions leading to a +337 pcm error in the sodium void worth. This last

) 5 . 6 .

 56 Compute axial leakage L i with eq. (23) Go back to point 2 until convergence 7. Compute effective cross sections Initialization ∀i,

Fig. 2 .

 2 Fig. 2. ASTRID internal fuel assembly: (a) axial layout, (b) fuel pins mesh, (c) sodium plenum mesh, (d) axial protection mesh

= 1 )

 1 without impacting the final results. This fact can be formally checked in the second line of

) 5 .

 5 Compute axial leakage L i with eq. (23) 6. Compute 2D fluxes ψ i with eq. (15)7. Compute effective cross sectionsIn this optimized version of the 2D/1D algorithm, 2D fluxes ψ i are used to weight cross sections in point 7. Equations[START_REF] Rahnema | Generalized energy condensation theory[END_REF] [START_REF] Douglass | Consistent generalized energy condensation theory[END_REF] and[START_REF] Vidal | APOLLO3 Homogenization Techniques for Transport Core Calculations -Application to the ASTRID CFV Core[END_REF] are used.

  Table II that shows the reactivity error when only overcritical layers are computed in iteration 1.

	Computed in	Iteration n o
	iteration 1	1	2	3
	All layers	+164 +69 +69
	Only over-critical	-6	+69 +69

  to obtain the standard flux volume weighting formula.

		Nominal	Voided	Void worth
	Reference	k ± 10 5 δk	k ± 10 5 δk	∆ρ Na (pcm)
	Monte Carlo 1.08183 ± 1 1.06399 ± 1	-1550 ± 2
	2D/1D	∆ρ (pcm)	∆ρ (pcm)	∆(∆ρ Na )
	φ lm	-58	-39	+19
	φ 00	+288	+625	+337

TABLE III

 III Results for fuel assembly calculation in nominal and voided situations. Reference eigenvalue / sodium void worth and reactivity discrepancy (towards reference) using 33 group homogenized cross sections generated with 2D/1D algorithm. φ lm (resp. φ 00 ) stands for flux moments (resp. scalar flux) averaged XS.

ACKNOWLEDGMENTS

The authors acknowledge the CEA of Cadarache for funding this work and AREVA and EDF for their long term support. They also thank the APOLLO3 R development team for their effort in developing the code. One of the authors (B.F.) would like to address many thanks to S.

Santandrea and E. Masiello for their great help while developing the methods presented in this paper.

value is comparable to typical SFR delayed neutron fractions β ≈ 360 pcm and thus far from being compatible with safety requirements. Therefore, the conclusion is that neutron transport physics in SFRs and even CFV-like designs can be accurately modeled with 33 group homogenized cross sections. Consistent information is however required for high angular order data else results strongly deteriorate.

Looking at the computational efficiency of the 2D/1D algorithm, a gain of a factor 5 to 10 was found on both calculation time and memory requirements compared to full 3D MOC simulations (with axial polynomial expansion [START_REF] Santandrea | Accelerated polynomial axial expansions for full 3D neutron transport MOC in the APOLLO3 code system as applied to the ASTRID fast breeder reactor[END_REF]). We however point out that results are not fully comparable because acceleration methods for the TDT-B solver are still under development, making it possible to expect a greater gain in the time needed to solve the 2D/1D equations. All calculations have been run with OMP parallelism on a 4 threads cluster node with Intel(R) Xeon(R) CPU E5-2630 v4 @ 2.20 GHz.

V. CONCLUSION

This paper presents an efficient algorithm to produce effective cross sections for SFR whole core transport calculations.

The 2D/1D method is used to approximate the 3D Boltzmann transport equation while preserving consistent representation of angular fluxes. The stability of the algorithm is numerically checked and the accuracy of the effective cross sections verified against reference Monte Carlo simulations.

The benchmark chosen for application is a radially reflected fuel assembly representative of the ASTRID CFV core. It is shown that very good accuracy is achieved with the 33 group homogenized 2D/1D cross sections, both in nominal and voided situations. It is also proven that high order angular information must absolutely by stored in the effective cross sections.

The algorithm presented in this paper is subject to two main sources of errors that should be commented. The first one is the flat flux approximation that is used to compute transverse leakage and homogenize 2D cross sections. If consequences of this assumption are likely to be limited in SFRs where neutrons have large mean free paths, it should be used with caution for other applications such as PWR calculations.

The second source of error is the fix-up used for negative sources. Because it replaces physical leakage by empirical ones and even neglects them in sub-critical layers, it breaks the neutronics balance in 2D calculations. Further work shall investigate quantitatively the impact of this approximation on effective cross sections to state whether it should be mitigated or not (for instance with a transverse leakage splitting technique).

Finally, undergoing work tends at confirming that the 2D/1D algorithm produces accurate effective cross sections for realistic industrial applications i.e. full core calculations with radial reflector, control rods. . . As for reactor cycle calculations, the algorithm provides a natural way to parametrize cross sections with realistic local burn-ups: the 3D fuel assembly power (thermohydraulics data) can be used to normalize 1D fluxes and run the isotopic depletion calculations for each axial evolving layer. Future work shall also investigate this possibility.


J Comte

C Guy

L Gosmain

S Parraud

Determining the porosity and water impregnation in irradiated graphite

Keywords: Graphite, nuclear waste, water uptake, radionuclides, porosity, UNGG

Impurities in nuclear graphite can become neutron-activated during operation, generating radionuclides and leading to modifications in the microstructure of the graphite due to fast neutrons. The long-term disposal of nuclear graphite waste requires collecting data on the behaviour of long-lived radionuclides in disposal conditions (water-saturated). The release of radionuclides in solution depends on a number of physicochemical processes such as water ingress into the graphite structure (radionuclides sites), as well as the solubilisation and transport of radionuclides in solution through the graphite pores. Analysis of the impregnation (impregnation speed, impregnation rate) of water in the porous graphite environment represents one of the main parameters that will greatly influence the physicochemical processes controlling the release of radionuclides in solution. For this reason, the impregnation of irradiated samples from the G2 and St Laurent A2 (SLA2) gas-cooled graphite-moderated reactors was studied, as was the distribution of the porosity. Results show that the geometric density of the samples decreases after irradiation, which is expressed as an increase in the total porosity or more precisely in the open porosity. This means that irradiated graphite is mainly a macroporous solid. Nuclear graphite becomes a more hydrophobic material with its open porosity filled by water in just a few days. This is due to the presence of functional groups on the pore surface.

Introduction

The decommissioning of graphite-moderated nuclear reactors will generate about 230,000 tonnes of graphite waste worldwide [START_REF] Banford | Strategic Options for the Management of Waste Irradiated Graphite[END_REF]. In France, the Commissariat à l'Énergie Atomique et aux Énergies Alternatives (CEA) and Electricité de France (EDF) operated 9 gas-cooled graphite-moderated reactors (UNGG) which have all been shut down for at least 20 years now. They were graphite-moderated and fuelled with natural metallic uranium. The decommissioning of these reactors will generate about 23,000 tonnes of irradiated-graphite (igraphite), representing around 81,000 m 3 of conditioned waste in cementitious packages [2]. Some of the impurities contained in the nuclear graphite became neutron-activated during operation, giving rise to numerous radionuclides such as 36 Cl, 14 C, 3 H, etc.

According to the French 2006 Act on the sustainable management of radioactive materials and waste passed by the French parliament [3], a dedicated shallow disposal concept was chosen as the reference solution for the long-term management of French graphite waste. In this concept of disposal site, the waste will be at term in contact with underground water. Thus, it is very important to know and understand underwater graphite behavior in terms of water uptake and hence of radionuclide release. In order to design the best long-term management option for irradiated -graphite, the producers (the CEA and EDF) carried out experimental studies in collaboration with the French national waste management agency (Andra) to determine the water uptake kinetic, the rate of water saturation and the radionuclide source term of irradiated graphite.

Natural graphite is a crystalline allotrope of carbon, made up exclusively of sp2 hybridised carbons. The graphite crystal belongs to the hexagonal system and consists in a compact stacking (AB stacks) of polycyclic aromatic layers (graphene layers). sp2 graphitic materials are most of the time considered as an hydrophobic materiel [START_REF] Kozbial | Are Graphitic Surfaces Hydrophobic?[END_REF].

Nuclear graphite used in UNGG or other graphite-moderated reactors (GCR) is a synthetic material very far from the crystalline structure of the natural graphite. For UNGG reactors, it was manufactured from petroleum coke mixed with coal-based binder pitch heated and extruded into bricks that were heat-treated up to about 2800/3000 °C for graphitization.

During this step, cleaning agents such as NaF or MgF 2 were used to obtain high-purity nuclear graphite grade [5-6-7]. As a result of the mixing of several carbon compounds, this material is structurally heterogeneous at a local scale (more or less anisotropic and not completely graphitized, with internal porosity….). Local heterogeneities in density may be caused either by blowholes or lack of impregnation of the open porosity during the manufacture of the graphite (local drop in density), or by the opposite by the presence of areas very poor in porosity, e.g. due to dense agglomerations of coke grains. Nuclear graphite is more a "carbon-carbon composites" relatively highly graphitised (blend of petroleum coke grains and less highly graphitisable coal-tar pitch) than a monocrystalline natural graphite.

The graphite blocks made with this material display good mechanical properties and chemical inertness. Some impurities remain at trace level in the nuclear graphite. Nuclear graphite contains small amounts of impurities like oxygen, hydrogen, metals and halogens, among them chlorine.

Irradiation of nuclear graphite but also temperature or radiochemical corrosion during reactor operating, have a major impact on the microstructure of the material at micrometric and nanometric scales [7-8 -9 -10]  The neutron flux induces displacement cascades of carbon atoms in the graphite matrix leading to structural defects in the material. Indeed, if the neutrons energy is sufficient (E> 0,1 MeV), the carbon atoms are ejected out of their initial position and can be for instance displaced between graphite planes , leading to creation of defects in the structure (interstitials, gaps). In the case of irradiation at low temperature, the consequences of the occurrence of these defects are the accumulation of stored energy (call Wigner energy).

 Irradiation temperature is also known to affect the grahite structure during irradiation.

A temperature between 200°C up to 500°C has been found to "heal" some defects, especially the small ones [11 12]. Defects can migrate and interstitial atoms could locally form new carbon layers.

 CO 2 radiolysis, mainly initiated through gamma irradiation, induces radiolytic corrosion processes at the gas/graphite interfaces with as consequence a growth of the pore size [6 -7 -13] and a decrease of the density.

CO 2(g) + hCO(g) + O*(g) O* + C(s) CO(g)
Which taken together is equivalent to the Boudouard reaction: CO 2(g) + C(s)2 CO(g) ΔH = + 172kJ mol -1

 To avoid or limit the radiolytic corrosion, corrosion inhibitors, such as methane, have been introduced into the coolant during operations, but their radiolytic decomposition leads to the build-up of undesired carbonaceous deposits. Recent analysis of irradiated graphite from two Magnox reactor cores highlighted the presence of a carbonaceous deposit on the exposed surfaces of the graphite bircks [START_REF] Heard | Evaluation of surface deposits on the channel wall of trepanned reactor core graphite samples[END_REF].

On a chemical point of view the modification of the C bonding induced by irradiation, radiolytic corrosion or by corrosion inhibitors on the surface lead to the creation of functional group [15 - 16 -17 -18]. More precisely, series of carbon-oxygen functional groups have been identified on the surface of irradiated graphite [16 -18 -19].

As a consequence the structure (porosity, microstructure and nanostructure) and the surface properties (presence of functional group) of the irradiated graphite are modified and the kinetic of water impregnation must be determined experimentally. That requires also an accurate knowledge of the irradiated graphite property. To assure a phenomena interpretation of water uptake the resulting porosity of the sample must be also determined. It must therefore be emphasized that the graphites used in the context of the UNGG reactor type are materials that must be more regarded as "carbon composite" (coke + binder) than as homogeneous materials.

The i-graphite of two different reactors -G2 and SLA2were investigated for this study.

These reactors were chosen due to the different characteristics of their graphite stack and their operating conditions, particularly the type of coke used for the graphite manufacturing, the neutron and thermal power levels, the CO 2 pressure and the thermal history of the graphite during reactor operation). Irradiated and unirradiated samples from both reactors were available for comparative studies.

Material and methods

Graphite samples

Graphite from the French UNGG G2 reactor G2 was the first gas-cooled graphite-moderated reactor designed to produce both electricity and plutonium. It was operated by the CEA from 1958 until 1980 at Marcoule in the south of France [20-21-22]. The graphite used in the G2 reactor was manufactured by Pechiney SA in

France. Purified and impregnated, this special coke graphite was used to build the active core of G2. This graphite can be considered as similar to "grade A" graphite previously used in the UK Magnox reactors. During operation, the graphite temperature ranged from 140°C to 380°C.

Within the framework of this study, the irradiated samples chosen were cored in the middle of the pile, as shown in Figure 1. Three of them were taken from the moderating part of the core (numbers 27, 32 & 42 in Figure 1). Unirradiated samples were also studied.

The main characteristics of the irradiated samples from G2 are shown in Table 1. The Saint-Laurent A2 (SLA2) reactor was commissioned in August 1971 and shut down in 1992. Its thermal power was 1700 MW. The gross electrical power was 530 MWe [START_REF] Bastien | Réacteurs à uranium naturel graphite-gaz (Uranium Natural Graphite-Gas reactors)[END_REF]. The graphite stack of the SLA2 reactor was cylindrical in shape with a vertical axis, a diameter of 15.73 metres and a height of 10.2 metres. The total mass of graphite was 2,440 tonnes. The graphite in the SLA2 stack (moderator or reflector) was made from LIMA coke, which underwent one impregnation and was purified using MgF 2 . All the graphite was manufactured by the Pechiney/SERS between May 1966 and November 1967. The operating temperature ranged between 240°C (top) and 470°C (bottom).

This study used twelve samples from two different fuel channels (F10M16-C20 and F7M15-C19) located at a radius of 5 m and 2.95 m respectively from the reactor centre with variable sampling heights (see Table 2 and Figure 2). Three samples per channel were chosen for structural characterisations and three for water impregnation tests.

FIGURE 2: SECTION THROUGH THE GRAPHITE STACK OF THE SAINT-LAURENT A2

The main characteristics of the G2/SLA2 reactor and graphite are given in Table 3

Origin and preparation of samples

Graphite from the French UNGG G2 reactor

Each sample was considered as homogeneous from a neutronic and thermal point of view.

Successive slices of about 1 to 3 cm were cut using a wire saw in a glove box (Figure 3). Some slices were re-cut to obtain smaller samples for porosity measurements. The dimensions of unirradiated and irradiated samples used for the water uptake tests are given in Table 5.

Other samples were used for porosity measurements. 4). In order to perform density and water uptake measurements, it was necessary to re-machine the samples. For density measurements, samples were cut into semi-circles of about 10 mm thick. For the six samples chosen for the impregnation tests, these semi-circles were used to measure porosity (helium pycnometry and mercury porosimetry). The rest of the core, i.e. about 40 mm long, was kept for continuous water impregnation tests. Figure 4 shows how the samples were divided up for the different characterisations. 6.

Materials and methods

Water uptake methodology

The study involved monitoring the change in the mass of samples immersed in water over time. Two techniques were used: continuous measurements for immersed samples suspended from a precision electronic balance (Test A), and intermittent measurements of the mass variations in immersed samples placed in experimental vessels (Test B).

In the first case, the sample was hung from a Sartorius balance beam and fully immersed in a container filled with ultrapure water. This configuration made it possible to continuously measure any changes in the apparent mass of the samples in water. As the water impregnates the sample, the buoyant force decreases and the apparent mass increases, which therefore makes it possible to determine the water mass having impregnated the sample. Measurements were taken over periods of about 10 to 40 days depending on the sample.

In order to determine the long-term impregnation kinetics, the samples were placed in reaction vessels filled with ultrapure water while making sure the S/V remained the same. The sample masses were measured regularly as follows: the sample was removed from the reaction vessel, gently wiped with paper towelling to remove the surface layer of water, weighed, and then put back in the reaction vessel.

Density measurements

Geometric densities were obtained by measuring the dimensions of the samples with a sliding calliper (uncertainty of 0.1 mm) by recording several readings per level, and weighing samples on a precision balance. The propagation of uncertainties, notably on average values, was calculated in accordance with standard engineering techniques.

Helium pycnometry

The non-destructive technique known as pycnometry by helium displacement is used to measure the internal volume of a porous material and to deduce its density from its mass. The measurement principle relies on the properties of helium, which behaves like an ideal gas. It is first pressurised within a reference chamber of a known volume before opening the pressurised chamber (sample) of a known volume. Pressure changes are associated with the porosity of the sample. Gas fills up the open porosities of the sample, while only the dense fraction and the sealed porosities of the sample remain inaccessible. The cells being used allow measurements on samples ranging from 1 to 10 g.

Mercury porosimetry

The destructive technique known as porosimetry by mercury intrusion is designed to determine the distribution of the pore sizes within a porous material (solid or powder). The measurement principle relies on the properties of a non-wetting liquid at ambient temperature, such as mercury. Mercury needs to be forced to penetrate the pores of a porous material. Any intrusion pressure, P, which is necessary for mercury to penetrate a pore with a diameter, d, is provided by the Laplace relation, as follows:

d P   cos 4  (1)
where P is the intrusion pressure of mercury (Pa)  is the surface tension of mercury (Pa/m)  is the wetting angle of mercury on pore surface (°)

d is the pore diameter (m)
The pressure applied to the mercury is increased by increments and the volumes of mercury penetrating the porous solid are recorded at each increment. The cumulative pore volume is then obtained as a function of the intrusion pressure or the diameter of the pores. The derivative gives the pore distribution of the solid. Knowing the density of the solid, the total quantity of mercury introduced into the solid can be used to determine the total porosity. This method of study has some drawbacks. In particular, the high pressure (about 400 MPa) applied to the mercury for it to penetrate the finest pores (a few nanometres in diameter) may alter the material studied. For this reason, the main advantage of this technique is its ability to determine the pore size distribution. The use of Laplace's law amounts to considering all the pores as cylinders, which is far removed from reality. Moreover, it must be assumed that the only sensitive factor is the diameter of the pore entrances. But significant hysteresis phenomena may conceal the actual pore structure. However, this technique is very useful for comparing different materials and provides a conventional representation of the pore structure.

During this study, a Micromeritics porosimeter (Model 9420 in a glove box,) was used to record measurements in a radioactive environment on inactive or irradiated samples.

Results and discussion

Part 1 -Density and porosity

Density

The results are shown in Table 7 for the G2 samples and in the Table 8 for the SLA2 samples.

As mentioned earlier, the G2 graphite cores were recut to obtain specific samples for water impregnation tests and specific samples were used for porosity measurements. The geometric density readings from prismatic and cylindrical lamellae proved to be consistent for all the irradiated samples (Table 7).

In case of the SLA samples, only one semi-circular disc was used for helium and mercury porosity measurements. Due to the lack of irradiated samples and the relatively small dimensions, it was not possible to perform reproducibility studies for the density and porosity measurements.

Sealed porosity

It can be seen that the irradiated graphite densities are mostly slightly lower than those of the unirradiated graphite. This decrease in density is associated with an increase in open porosity, and is certainly related to the wear of graphite under irradiation (radiolytic corrosion by the CO 2 coolant).

There are some exceptions: sample No. 42 from G2, whose density can be considered as identical as before irradiation, was the first moderator core sample taken at the limit of the reflector. The porosity values are almost identical to those determined for the inactive brick. It was not as irradiated and was exposed to a lower temperature than the two others. This leads us to believe that radiolytic corrosion may be considered as negligible in this area (at the boundary between the reflector and the moderator).

Two SLA2 samples (SLA2-58 and 60) revealed an equivalent or higher density than those determined for unirradiated graphite, together with a very low open porosity. These samples must have come from a batch manufactured with a higher level of impregnation.

In the other case of the G2 or SLA moderator samples, it could be seen that the increase in the total porosity was linked to an increase in the open porosity. More precisely, it appears that the open porosity level generally correlates with the neutron fluence received by the samples, and therefore with the position in the graphite stack. The level of closed porosity remains constant or slightly higher than the reference in the case of SLA2 samples, whereas the sealed porosity decreased in the case of G2. It is more complex to assess the variation in closed porosity according to the position of the samples in the core. In light of these results, the decrease in sample density definitely results from the radiolytic corrosion of the graphite.

Pore size distribution

The variation in pore volume according to the diameter of the pores and the pore distribution of the samples are shown in Figure 5 for unirradiated samples, in Figure 6 for the G2 sample, and in Figure 7 for SLA2 samplesErreur ! Source du renvoi introuvable.. It can be seen that the pore size distribution for unirradiated samples differs between the two types of graphite. For G2 graphite (special coke graphite), around 75% of the total porous volume of the sample lies within a pore diameter range between 1 and 30 µm with a maximum distribution oscillating towards 3.5 µm.

In the case of the SLA2 samples manufactured with LIMA coke, it appears that unirradiated graphite has two types of distribution: around 35% of the pores in the mesoporous domain (pore diameter between 1 and 300 nm), and 40% of the total porosity in the macroporous zone between 1 to 4 μm with a peak at around 2.2 µm. These results are in accordance with some characterisations performed on other types of nuclear graphite [START_REF] Li | The Influence of Pores on Irradiation Property of Selected Nuclear Graphites[END_REF]. A general increase in the macroporosity could be seen in the G2 graphite after irradiation (pore diameters ranging between 1 and 30 µm), including a wider pore size distribution mostly for the G2-27 (68% porosity) and G2-32 (72% porosity) samples with a modification in the average pore radius in the distribution. The average radius increased from 3.5 µm on inactive samples to 4.3 µm for sample G2-42, then to 5.4 µm for sample G2-32, before reaching 5.7 µm for sample G2-27.

Overall, it appears that the fraction of pores in the mesoporosity domain (pore diameter between 1 and 300 nm) is not affected in the case of SLA2 irradiated samples. The intensity of the peak in porosity at about 1-4 μm decreases. Finally, there is a third porosity zone around 10 μm, but with greater scatter (varying from 5 to 20 μm) according to the samples. Irradiated samples from two reactors and with different neutron flux or temperatures were characterised, as were unirradiated samples. Generally speaking, irradiation leads to a relatively minor change in the characteristics of the graphite, the average wear being less than 5%, although it reaches close to 8% locally. Overall, a decrease in the density could be observed, associated with an increase in open porosity mainly in the macroporous domain.

The results are not significant enough to distinguish any general variations in closed porosity.

These changes in the characteristics are associated with the radiolytic corrosion of graphite, corrosion which is all the greater when the samples are subjected to high fluences and irradiation temperatures.

Part 2 -Water up-take

Unirradiated samples

The water uptake kinetics for unirradiated samples are presented in Figure 8. It can be seen that the water uptake kinetics are slower for the G2 samples than for the SLA2 samples.

Moreover, only less than half of the total open porosity, even after more than 700 days of impregnation, was filled for the G2 samples. Unirradiated SLA2 graphite presents a slightly different behaviour. The water impregnation rate is faster but, at the end, water uptake seems to be limited to around 70% of the total porosity. The water impregnation kinetics for the two types of graphite were relatively high during the first phase lasting about 150 days and then they slowed down dramatically.

A connection can be established between these two steps of the impregnation process and the pore distribution of nuclear graphite using mercury intrusion porometry. The first rapid-filling step relates to the filling of the macroporosity that represents 65% of the total open porosity in the case of the SLA2 graphite samples (35% of the pores are in the mesoporosity domain pore diameter between 1 and 300 nm). For G2 graphite, the water does not fill all the macroporosity of the sample. This difference cannot only be attributed to the diameter of the pores but also to the hydrophobic phenomena. To check this assumption, some complementary experiments were performed in an ethanol-water mixture. Adding ethanol (10 wt%) improves impregnation: the kinetics and saturation rate increased significantly. Adding ethanol to water helps lower the surface tension (Table 9) of the mixture and therefore promotes its impregnation in a hydrophobic solid. In the ethanol-water mixture, water uptake in the G2 samples is faster as illustrated in Figure 9 versus the square root of time, and the saturation of the open porosity increases from 40% to 60%. To conclude, the water uptake in unirradiated graphite not only depends on the pore size distribution, but also on the initial surface state.

These differences between nuclear graphite from different origins are reduced when considering irradiated samples. Irradiated graphite is a more hydrophilic material than the unirradiated one. It is evidenced in Figure 10 that the porosity is rapidly filled in the first 30 to 50 days and then this slows down to reach 75% to around 100% of the open porosity.

At the same time, the curve in Figure 10 shows that the difference in size between the G2 samples does not seem to modify the total rate of open porosity, but only the water impregnation kinetics.

These results show that irradiation increases the water impregnation kinetics and the saturation rate of the graphite. The original behaviour of unirradiated graphite is completely erased by the irradiation time in the reactor. Two reasons can explain this difference:

 First, irradiation modifies the macroporosity of the different types of graphite, particularly by widening the mean pore diameter as shown in the mercury intrusion porosimetry spectrums, which may facilitate water impregnation. Irradiated graphite was mainly a macroporous solid (about 80% of the total porosity with an apparent pore diameter ranging between 1 and 30 µm).

 Second, it is known that irradiation and radiolytic corrosion result in breaking the C-C bonds in the graphite crystallites and/or grain boundaries to form more hydrophilic C-H or C-O bonds at the surface.

These two phenomena lead to faster and higher impregnation.

Conclusion

In the framework of the studies on the disposal of graphite waste from UNGG reactors, it is necessary to obtain reliable data on the behaviour of long-lived radionuclides when graphite will be water-saturated in disposal conditions. Water impregnation in irradiated graphite is the first step of a process that leads to the release of radionuclides contained in the graphite. If this step is slow and low, it can control both the kinetics and the release rate of radionuclides in solution. Impregnation tests were therefore performed on unirradiated and irradiated graphite samples from the G2 and St Laurent A2 reactors which had been previously characterised by geometric density measurements, He pycnometry and Hg porosimetry.

Characterisation of the structural (or microstructural) properties and the water uptake experiments of unirradiated graphite made it possible to determine that they do not have the same pore size distribution and most probably have a different surface effect depending on their origin. The SLA2 graphite was seen to have two types of porous distribution, i.e. around 35% of pores in the mesoporous domain (pore diameter between 1 and 300 nm) and 40% of the total porosity in the macroporous zone (pore diameter between 1-5 μm). G2 graphite presents a quasi-single macroporous distribution. Moreover, SLA2 graphite appears to be more hydrophilic than G2 graphite.

In general, irradiation induces a small modification in the characteristics of graphite. In both cases, it appears that the geometric density of the samples decreases under irradiation, which is expressed as an increase in total porosity, more precisely in the open porosity. The pore size in the macroporous domain increases (total fraction and average diameter), whereas the mesoporous fraction does not seem to be affected by irradiation. This variation appears to be correlated with the neutron flux. This is a consequence of radiolytic corrosion of the graphite under CO 2 irradiation. Results from water uptake experiments show that the graphite becomes more hydrophilic after irradiation. The original behaviour of unirradiated graphite seems to be completely erased by their time in the reactor. This is probably due to the combined effects of temperature and irradiation, which lead to radiolytic corrosion of graphite and the presence of functional groups on the surface, which seems to modify the surface tension. The impregnation kinetics are faster and the saturation rate is close to 100% with regard to open porosity.

Considering the leaching behaviour of radionuclides in nuclear graphite under disposal conditions, these results show that the open porosity can be considered as being filled with water in a hundred days. Water penetration can be considered as a limiting factor for the most mobile radionuclides (leaching under <100 days), but not for the non-mobile radionuclides.

Moreover, the porosity of irradiated graphite in storage conditions will be saturated by the water site which is alkaline water. A possible chemical degradation in the surface compounds by this aggressive media is to be considered.

FIGURE 4 :

 4 FIGURE 4: SLA2 GRAPHITE SAMPLES AND THEIR DISTRIBUTION FOR THE DIFFERENT ANALYSES

FIGURE 5 :

 5 FIGURE 5: COMPARISON OF PORE DISTRIBUTIONS FOR UNIRRADIATED GRAPHITE SAMPLES FROM THE SLA2 REACTOR MODERATOR (LIMA COKE) AND THE G2 REACTOR MODERATOR (SPECIAL COKE)

FIGURE 6 :

 6 FIGURE 6: PORE DISTRIBUTIONS FOR IRRADIATED AND UNIRRADIATED GRAPHITE SAMPLES FROM THE G2 REACTOR MODERATOR (SPECIAL COKE)

FIGURE 7 :

 7 FIGURE 7: PORE DISTRIBUTIONS FOR IRRADIATED AND UNIRRADIATED GRAPHITE SAMPLES FROM THE SLA2 REACTOR MODERATOR (LIMA COKE)

FIGURE 9 :

 9 FIGURE 9: COMPARISON OF WATER AND WATER-ETHANOL IMPREGNATION RATES FOR UNIRRADIATED STACK GRAPHITE FROM THE G2/SLA2 REACTORS

TABLE 1 :

 1 IRRADIATED SAMPLES FROM THE G2 MODERATOR

	Sample No.	Sampling height in the reactor (m)	Thermal history (°C)	Fast fluence E>0.1 MeV (n.cm -2 )	Initial mass (g)	Dimensions (mm)
	G2-27	13.60-13.80	327	4.7.10 21	942	H=160 ø=63
	G2-32	14.60-14.80	320	4.7.10 21	863	H=150 ø=63
	G2-42	16.60-16.80	309	4.0.10 21	845	H=150 ø=63

TABLE 2 :

 2 IRRADIATED SAMPLES FROM THE SLA2 MODERATOR

	172
	Sample No.	Channel	Sampling height in the reactor (mm)	Thermal history (°C)	Fast fluence E>0.1 MeV (n.cm -2 )	Initial mass (g)	Approximate dimensions (mm)
	SLA2-43	F10M16 -C20	1680	430	~ 2.0.10 21		H=50 ø=19
	SLA2-44	F10M16 -C20	2070	435	~ 2.3.10 21		H=50 ø=19
	SLA2-53	F10M16 -C20	6500	330	~ 3.0.10 21		H=50 ø=19
	SLA2-55	F10M16 -C20	7280	310	~ 2.8.10 21		H=50 ø=19
	SLA2-58	F10M16 -C20	8660	270	~ 1.9.10 21		H=50 ø=19
	SLA2-60	F10M16 -C20	9260	250	~ 1.4.10 21		H=50 ø=19
	SLA2-122	F7M15 -C19	1080	440	~ 1.5.10 21		H=50 ø=19
	SLA2-124	F7M15 -C19	2070	455	~ 2.3.10 21		H=50 ø=19
	SLA2-129	F7M15 -C19	4480	400	~ 3.3.10 21		H=50 ø=19
	SLA2-135	F7M15 -C19	7280	310	~ 2.8.10 21		H=50 ø=19
	SLA2-138	F7M15 -C19	8660	270	~ 1.9.10 21		H=50 ø=19
	SLA2-139	F7M15 -C19	9260	250	~ 1.4.10 21		H=50 ø=19
	173

Table 4 TABLE 3 :

 43 and OPERATING CONDITIONS OF THE G2 AND SLA2 REACTORS

	Reactor	G2	SLA2
	First criticality	July 1958	June 1971
	Shutdown	February 1980	May 1992
	Thermal power (MW)	260	1700
	CO 2 pressure (MPa)	1.5	28.5
	Mass of graphite stack (tonnes)	1,500	2,200
	Graphite temperature during operation	140-380°C	240-470°C

TABLE 4 :

 4 CHARACTERISTICS OF THE G2 AND SLA2 REACTOR GRAPHITE (BEFORE

		IRRADIATION)
	Reactor	G2 Moderator	SLA2
	Coke	Special A grade coke NaF cleaned	LIMA coke MgF 2 cleaned
	Graphite density	~1.71	~1.68
	Capture section (mbarn)	3.95	3.75
	Thermal expansion coefficient (25-525°C) (//)	1.25.10 -6 K -1	2.4.10 -6 K -1
	Anisotropy ratio (┴)/(//))	2.3	1.6
	Compression (//)	27 MPa	42 MPA
	Ash content (ppm)	106 to 125	98

TABLE 5 :

 5 GEOMETRICAL CHARACTERISTICS OF THE SAMPLES FROM G2 FOR WATER

	UPTAKE TESTS

10 G2-32-sample FIGURE 3: CUTS MADE ON THE CORE OF THE G2 GRAPHITE MODERATOR  Graphite from the French UNGG SLA2 reactor Initially irradiated samples were in the form of cylindrical cores (approximate dimensions: diameter 19 mm / height 50 mm) with a bevelled end corresponding to the outer surface of the brick (Figure

TABLE 6 :

 6 GEOMETRICAL CHARACTERISTICS OF THE SAMPLES FROM SLA2

		Water impregnation
						Hg & He porosimetry
		19 mm
			37 mm minimum	 3-5 mm
			Reference			Mass	Average dimensions	Surface	Volume
	Graphite	sample	form	Water uptake test	(g)	Diam.*L. (mm*mm)	cm 2	cm 3
	Unirradiated graphite	Z1 Z2 XY1	cylinder cylinder cylinder	A A A		4.0085 7.7454 7.5475	12*21 12*41 12*39	10.1±0.3 2.35±0.03 17.7±0.4 4.63±0.02 17.1±0.4 4.44±0.02
		SLA2-44 cylinder	A		21.9988	19*48	34±1	13.6±0.4
		SLA2-55 cylinder	A		17.9422	19*39	29±1	11.1±0.3
	Irradiated	SLA2-58 cylinder	A		22.1746	19*47	33.5±0.9 13.2±0.3
	graphite	SLA2-124 cylinder	A		22.0483	19*50	35.7±1.1 14.2±0.4
		SLA2-135 cylinder	A		23.1589	19*51	36.1±1.9 14.5±0.8
		SLA2-138 cylinder	A		22.0355	19*50	35.4±2.7 14.1±1.1

TABLE 7 :

 7 DENSITY AND POROSITY OF THE SAMPLES FROM G2

		Average density ( G )	 S Helium	Total porosity Ref d=2.266	Open porosity (1- G/  S) )	Sealed porosity
	UNIRRADIATED
	Coke special	1.69 ± 0.06	2.15 ± 0.05	25.9 ± 0.9	21.7 ± 1.2	4.9 ± 1.4
	(average n=8)
	G2-27-1	1.67±0.05	2.17±0.09	26.3±0.6	22.9±1.1	3.4±1.2
	G2-27-3	1.62±0.04	2.18±0.05	28.5±0.4	25.7±1.3	2.8±0.8
	G2-27-7	1.65±0.05	2.27±0.09	27.2±0.7	27.3±1.3	-0.1±1.4
	G2-27-8	1.68±0.05	2.22±0.09	25.9±0.6	24.3±1.1	1.7±1.3
	G2-27 MEAN	1.66±0.06	2.21±0.09	27.0±0.6	25.0±1.1	1.9±1.3
	G2-32-1	1.67±0.05	2.16±0.09	26.5±0.7	22.91±1.1	3.5±1.3
	G2-32-3	1.62±0.05	2.19±0.05	28.5±0.7	26.0±1.0	2.5±0.8
	G2-32-6	1.65±0.05	2.18±0.09	27.1±0.7	24.3±1.1	2.8±1.3
	G2-32-7	1.67±0.05	2.17±0.09	26.4±0.7	23.0±1.1	3.4±1.3
	G2-32 MEAN	1.65±0.05	2.18±0.09	27.1±0.7	24.0±.1.1	3.1±1.5
	G2-42-1	1.71±0.05	2.11±0.09	24.5±0.7	19.1±0.9	5.4±1.2
	G2-42-3	1.71±0.05	2.18±0.09	24.5±0.6	21.6±0.9	3.0±0.6
	G2-42-6	1.69±0.05	2.17±0.09	24.3±0.7	20.8±1.0	4.9±1.1
	G2-42-7	1.70±0.05	2.14±0.09	24.8±0.7	20.5±1.0	5.8±1.1
	G2-42 MEAN	1.70±0.06	2.15±0.09	24.5±0.7	20.5±1.0	4.8±1.1

TABLE 8 :

 8 DENSITY AND POROSITY OF THE SAMPLES FROM SLA2

	AVERAGE DENSITY ( G )	 S Helium	Total porosity Ref d=2.266	Open porosity (1- G/  S)

TABLE 9 :

 9

	Solution	Surface tension (mN/m)
	Water	72.8
	Ethanol	22
	Water+10% wt Ethanol	52

SURFACE TENSION AT 20°C FIGURE 8 : COMPARISON OF WATER AND WATER-ETHANOL IMPREGNATION RATES FOR UNIRRADIATED STACK GRAPHITE FROM THE G2/SLA2 REACTORS Considering the uncertainty on the measurement of the sample volumes, open porosities and mass gains, a relative uncertainty of 10% was attributed

Acknowledgements

This work was carried out under the CARBOWASTE Programme: Treatment and Disposal of Irradiated Graphite and Other Carbonaceous Waste (Grant Agreement number FP7-211333), with partial funding from ANDRA, EDF and the CEA. The authors would also like to thank those in charge of G2 or EDF reactor operations and dismantling for having provided these samples and Mr Laurent Petit, EDF, for his assistance and follow-up of these studies.

Hexagonal cell F10M16

Hexagonal cell F7M15

Irradiated samples

The results obtained on irradiated graphite are presented in Figure 10