O L Krutkin

A B Altukhov
email: a.altuhov@mail.ioffe.ru

A D Gurchenko

E Z Gusakov

M A Irzak

L Esipov

A V Sidoro

L Chôné

S P Kiviniemi

Leerink

P Niskala

C Lechte

S Heuraux

G Zadvitsk

Validation of full-f global gyrokinetic modeling results against the FT-2 tokamak Doppler reflectometry data using synthetic diagnostics

Keywords: gyrokinetic, Doppler reflectometry, synthetic diagnostic, radial correlation, phase modulation

Two versions of the X-mode Doppler reflectometry (DR) synthetic diagnostics are developed within the framework of the ELMFIRE global gyrokinetic modeling of the FT-2 tokamak ohmic discharge. In the 'fast' version the DR signal is computed in the linear theory approximation using the reciprocity theorem, utilizing the probing wave field pattern provided by computation and taking into account the 2D plasma inhomogeneity effects; whereas the alternative 'slow' version DR synthetic diagnostic is based on the full-wave code IPF-FD3D describing the probing and scattered wave propagation in turbulent plasma. The DR signal frequency spectra and the dependence of their frequency shift, width and shape on the probing antenna position are computed and shown to be similar to those measured in the high-field side probing DR experiment at the FT-2 tokamak. The geodesic acoustic mode characteristics provided by the measurements and by the synthetic DR are close within a 12% accuracy. However, a substantial difference was found in the decay of the DR signal cross-correlation functions with growing frequency shift in the probing wave channels. The quick decrease in the radial correlation DR coherence observed in the experiment and full-wave synthetic diagnostic, compared to the fast synthetic DR, is attributed to the nonlinear effect of the probing wave phase modulation by the turbulence in the former two cases. The variation in the DR signal at a growing incidence angle in the experiment is also shown to be slower than predicted by both of the synthetic diagnostics, presumably due to underestimation of the probing wave phase modulation and consequent nonlinear saturation of the DR signal at lower incidence angles in modeling.

Introduction

Anomalous energy transport remains one of the main unresolved enigmas on the way to the realization of magnetic confinement fusion using the tokamak concept. According to the present day understanding [1,2] anomalous transport is determined by multi-scale drift-wave turbulence and is strongly dependent on the nonlinear interaction of its components, such as large-scale mean E × B flows, meso-scale zonal flows and fine-scale micro-turbulence excited due to specific profiles of the plasma parameters. This interaction controlling anomalous transport has been the focus of theoretical and experimental research in magnetically confined plasmas for a long time. The massively parallelized particle-in-cell and semi-Lagrangian simulations of the gyrokinetic (GK) distribution function and the electric field [3,4] provide an efficient theoretical tool for studying nonlinear turbulent plasma dynamics, which, however, needs validation and comprehensive benchmarking against the experimental data. In recent papers [5,6] such simulations of the electron and ion distribution function from the first principles were performed for the small research limiter tokamak FT-2. The complex interplay of the multi-scale turbulence components was studied in simulations, and successful quantitative comparison was presented with measurements carried out at the low magnetic field side of the machine with O-mode Doppler reflectometry (DR) and enhanced microwave scattering in the upper hybrid resonance [5,6].

In the present paper the results of the global GK particlein-cell simulations will be quantitatively compared to the X-mode DR experimental data obtained at the high magn etic field side of the FT-2 tokamak (these data characterize the tokamak turbulent dynamics and transport phenomena at different spatial and temporal scales). DR is the microwave diagnostic widely used nowadays in toroidal fusion experiments to obtain information on the multi-scale anomalous transport phenomena-namely, on the plasma rotation (both the mean value and its oscillations) and the turbulence poloidal wavenumber spectrum [7][8][9][10][11]. Recently, the radial correlation (RC) modification of this technique [12], providing information on the turbulence radial structure, was justified [13,14].

Unfortunately, the interpretation of DR data is complicated due to the contribution of the poorly localized smallangle scattering along the wave trajectory, which can lead to overestimation of the turbulence RC length in the linear scattering regime [13,14], and to its underestimation in the strongly nonlinear regime [15]. In the latter case the turbulence wavenumber spectrum measurements are questionable, and only the plasma turbulence mean velocity could be determined (however, with poor spatial resolution) [16]. Moreover, in the X-mode DR the turning point cannot be determined unambiguously due to 2D propagation effects, which additionally complicate the experimental data interpretation and their comparison with the theory predictions.

Development of a synthetic diagnostic that allows direct signal computation based on the results of GK modeling could be helpful in this case, both for interpretation of the experimental results and for the code benchmarking [17,18]. In this paper two versions of the X-mode DR synthetic diagnostics are compared within the framework of the ELMFIRE global GK modeling of the FT-2 tokamak ohmic discharge [6]. The X-mode DR signal is computed both in the linear theory approximation using the reciprocity theorem, utilizing the probing wave field pattern provided by computation and taking into account the 2D plasma inhomogeneity effects [19], and by the full-wave code IPF-FD3D describing the probing and scattered wave propagation in turbulent plasma [20]. The temporal behavior of the DR signal frequency spectra and the dependence of its amplitude, frequency shift and shape on the probing antenna position are computed and compared to those measured in the experiment at the FT-2 tokamak. In the case of multi-frequency probing the RCDR cross-correlation function (CCF) is also determined by both of the synthetic diagnostics and is compared to that obtained in the experiment.

The experimental approach

The experiment was performed at the FT-2 tokamak (the major radius R = 55 cm, the limiter radius a = 7.9 cm) in the hydrogen ohmic discharge (with plasma current I p = 19 kA, central density n e (0) = 4 × 10 13 cm -3 and electron temperature T e = 470 eV). The discharge is similar to that utilized for successful comprehensive benchmarking of the ELMFIRE GK code in [5,6], however, the toroidal magnetic field at the discharge axis is slightly smaller (B t (0) = 1.7 T instead of 2.1 T). The actual measured electron density and temperature profiles for this discharge used in the GK modeling, as well as the ion temperature profile, are shown in figure 1. The vertically movable (by ±2 cm) X-mode double antenna set (shown in figure 2), installed at the high magnetic field side in the vicinity of the equatorial plane, allowed plasma probing at a variable incidence angle with frequencies in the range f i = (50-75) GHz. It was used both in standard DR measurements utilizing the quadrature scheme and in the RCDR technique based on the scheme described in [21]. In the former the DR spectrum was obtained by averaging 155 random DR power spectra obtained using Fourier transform of 12.8 microsecond samples of the DR signal. In the latter the reference channel generator was tuned to the master frequency f 0 = 70 GHz determining the measurement position in the vicinity of r = 5 cm, whereas another generator (producing the probing signal at a slave frequency varied in the interval ±1.2 GHz around the master frequency by 100 MHz steps) was used in the second (signal) channel to determine the turbulence twopoint CCF.

The backscattering (BS) spectra P s (f ) (f = f sf i , where f i and f s are the probing and scattering frequencies, accordingly) measured by the X-mode DR diagnostics at different vertical antenna displacements were presented in [19]. It was shown there that the BS power decreases with the growing antenna vertical shift corresponding to the increase of the probing wave incidence angle, whereas the spectra frequency shifts and width increase.

The GK computation approach and the synthetic diagnostics

The simulations are performed with the global electrostatic particle-in-cell code ELMFIRE [3] to model the evolution of the full distribution function of drift kinetic electrons and GK protons and O 6+ ions. The simulation geometry has a circular cross-section and covers the whole radial range. The time step is set to δt = 30 ns and the spatial grid to 110 × 8 cells in the radial (r) and toroidal (ϕ) directions, while the number of cells in the poloidal (θ) direction is set to increase with r up to 950 at r = 7 cm. The corresponding computational grid cell sizes at the minor radius r = 5 cm are as follows: 0.07 cm radial, and 0.07 cm poloidal. They are sufficient for an adequate description of fluctuations contributing to the DR signal in the experiment (the poloidal wavenumber is 0.4 cm -1 ) and correspond to the measured RC length (radial wavelength is 0.7 cm). In the present version of the code [22] the parallel derivatives are calculated along the field lines (not as a combination of poloidal and toroidal derivatives at a given point). In practice this requires interpolation in the poloidal plane as the end points used in the derivative calculation are not grid nodes. As Hariri and Ottaviani [23] point out, with such a method, the coarse grid in the toroidal direction can be compensated by the fine grid in the poloidal direction carrying the information about the fine structure of the turbulence. The toroidal mode number is only the projection of parallel and poloidal components. However, modeling of smaller-scale turbulence, in particular, belonging to the electron temperature gradient (ETG) mode, is not possible with the ELMFIRE GK code. Approximately 2400 electrons and ions are simulated per grid cell. A momentum and energy conserving binary collision operator [24] is applied, allowing the inclusion of neoclassical physics in the simulations.

The particles are initialized according to a prescribed temper ature and density profiles taken close to the experimentally measured ones (see figure 1). The profiles are allowed to develop self-consistently in time, as turbulence develops and heat sources and sinks, namely ohmic heating and radiation losses, are applied. The ohmic heating is induced by a spatially homogeneous loop voltage [25,26], while electrons are cooled using a Monte Carlo model according to an impurity radiation loss distribution measured by the bolometry diag nostics. The initial oxygen O 6+ impurity density is chosen to provide the experimental value of Z eff = 2.2, as determined using ASTRA modeling and experimental loop voltage measurements, with equal temperatures for the impurities and main ions. Neumann and Dirichlet boundary conditions are used at the boundary of the simulation domain for the electrostatic potential. Particles passing the boundary are returned to the simulation domain as an electron-ion pair according to a probability distribution proportional to the measured radial profile of neutral hydrogen density and at the wall temperature.

The fast synthetic X-mode DR diagnostic developed for this code [19] is based on the reciprocity theorem of electrodynamics [27], relating the high-frequency current in the plasma volume and the signal radiated by it and received by the antenna. The probing and BS wave propagation in the plasma not perturbed by turbulence is described in the fast synthetic diagnostic by the full-wave code WaveTOP2D [28]. The code solves complete Maxwell equations in the toroidal geometry accounting for the magnetic field shearing and presence of metallic walls, however, it suppresses reflections from the tokamak walls by specially introduced damping in the surface layer. The antenna beam is assumed Gaussian, as it was measured and computed at the air until the -10 dB power level. Turbulent density fluctuations computed in the eight toroidal cross-sections were used to obtain the DR signal with the total duration 0.96 ms. This realization was divided into 0.012 ms samples in which the random Fourier power spectra were calculated and then used to obtain the average DR power spectrum. It should be mentioned, however, that the approach based on the reciprocity theorem neglects the perturbations of the probing wave propagation caused by the density fluctuations; in particular, multiple small-angle scattering or its strong phase modulation. In this sense it is equivalent to the so-called Born approximation linear in the fluctuation amplitude. The main merit of this approach is the possibility of fast computations, whereas the main drawback is related to neglecting the nonlinear effects.

The alternative 'slow' DR synthetic diagnostic is free from this drawback. It is based on the full-wave code IPF-FD3D describing the probing and scattered wave propagation in turbulent plasma [20]. It accounts for the probing wave perturbations introduced by density fluctuations and, therefore, is not limited to the linear approximation. However, it solves the Maxwell equations with perfectly matched layer boundary conditions for a single polarization mode in the case of the propagation perpendicular to the magnetic field, the poloidal component of which is neglected as well as the toroidicity. It should be mentioned that multiple reflections between the plasma and the tokamak wall or microwave antenna, which could be important for the DR signal formation at small probing wave incidence angles, are not taken into account by either of the synthetic diagnostics.

It should be noted that the full-wave synthetic diagnostic simulations are very time demanding. With the full-wave code the computation of the synthetic reflectometry signal for a specific set of FT-2 tokamak experimental conditions requires around four hundred hours. This time scales with linear size as (a/λ) 3 , where a is the plasma size and λ is the characteristic wavelength. The analogous computation by the fast synthetic diagnostic requires approximately four hours, but scaling with size is less favorable in this case-(a/λ) 4 . In both cases computations are performed on a single Intel Xeon 2.4 GHz CPU. Although full-wave synthetic diagnostics offer the possibility of parallelization and multi-frequency computation, they were not enough to overcome the drastic difference in computation time for the FT-2 experimental conditions.

DR frequency spectra comparison

The DR spectra obtained in the experiment and provided by the synthetic diagnostics are compared in figure 3 for vertical antenna shifts of different signs and values.

As shown in the figures, the spectra provided by both of the synthetic diagnostics fit reasonably well to the experimental ones for all the antenna positions. Only at the maximal antenna shifts do both of the synthetic diagnostics fail to describe the wing of the experimental spectrum opposite to the spectrum frequency shift direction. Nevertheless, not only are the shifts of both spectra close in all cases, but their shapes are also similar. This impression is confirmed by figures 4 and 5, where the mean frequency shift f D = ´f • P s ( f ) df / ´Ps ( f ) df and the mean frequency

width ∆f = î ´( f -f D ) 2 P s ( f ) df / ´Ps ( f ) df ó 1/2
of the experimental and synthetic spectra obtained as the first and the second moment of the spectrum, respectively, are plotted as functions of the antenna vertical displacement and the turbulence poloidal wavenumber. The relation between the antenna displacement and the fluctuation poloidal wavenumber was obtained using the X-mode ray tracing (as the doubled value of the probing wave poloidal wavenumber at the trajectory minimal radial coordinate). The spectral shift data provided by the experiment and GK modeling agree surprisingly well. A meaningful difference is only observed at the largest antenna displacement y a = ±2.0 cm. This difference is probably due to a contribution of experimental spectrum wings becoming especially heavy and asymmetric at high antenna displacement and, correspondingly, the fluctuation poloidal wavenumber. As was mentioned, the wing is opposite to the mean frequency shift of the spectrum and, presumably, is caused by small-scale drift modes excited in the plasma at a low level. These modes rotating in the direction of the ion diamagnetic drift might be poorly described by the GK code, since the smallest scales can be close to the limit of grid resolution and are also most sensitive to the particle noise. Furthermore, one should also note that in [21] for similar parameters we showed clear frequency mismatch, even with different signs for high poloidal wave numbers close to the edge when comparing the filtered ELMFIRE signal to the linear GS2 analysis. Since the disagreement decreases when increasing the plasma current it was assumed to be caused by wide orbit effects, which were not included in the GS2 analysis. However, since the mode spectrum is very sensitive to the small changes in gradient lengths such disagreement can also be a result of the small changes in profiles during the full-f simulation. Another possible explanation of this difference could be associated with a non-Gaussian form of the experimental antenna beam at the -10 dB level. Figure 4 allows us to obtain estimations for the fluctuation poloidal velocity in the experiment and in the computation. The estimation of the mean fluctuation poloidal velocity in the experiment is given by v θ = πf D /k θ ≈ 2.1 ± 0.2 km s -1 , whereas the synthetic DR computations result in v θ ≈ 2.0 ± 0.1 km s -1 for the fast diagnostics and v θ ≈ 2.4 ± 0.1 km s -1 for the full-wave one. The obtained agreement between the measured and computed DR spectra, which are determined by the plasma flows, provides evidence in favor of correct modeling of plasma poloidal rotation and a correct reproduction of the radial electric field behavior in the FT-2 tokamak by the ELMFIRE code.

As shown in figure 5, the mean width of the synthetic spectra is systematically smaller than in the experiment. However, the difference is in the range of only 20%-30%. It could be attributed to the contribution of the spectral wings, which are higher in the experiment.

The large backscattering spectrum width observed in figure 5 is partly associated with strong plasma potential oscillations caused by the geodesic acoustic mode (GAM) excited in the FT-2 plasma [5,6,29]. These oscillations could be extracted from the DR signal using information on the temporal behavior of the frequency shift of its spectrum [7]. By utilizing 5 µs samples of the experimental and synthetic DR signal the following dynamics of the plasma poloidal velocity shown in figure 6 were obtained. The oscillations in the v θ (t) behavior are caused by the GAM. The velocity oscillation amplitudes provided by the experimental measurement, fast and full-wave synthetic diagnostics are close (1.7 ± 0.3 km s -1 ; 1.6 ± 0.3 km s -1 and 1.4 ± 0.2 km s -1 , respectively). The amplitude of the plasma poloidal velocity oscillations determined directly from the GK code output using the magnetic surface averaged radial electric field is also close to these values (1.3 ± 0.5 km s -1 ). The standard deviations of the velocity fluctuations determined using DR data are also similar: the values 1.1 km s -1 , 0.9 km s -1 , 1.0 km s -1 and 0.8 km s -1 were obtained in the experiment, given by the two synthetic diagnostics and directly provided by the GK code. The GAM is represented by a line in the Fourier spectrum of velocity shown in figure 7. As shown in the figure, the experimental value of the GAM frequency is 40 kHz, and the computational value provided by the fast synthetic diagnostic is slightly higher at 45 kHz; whereas the value given by the full-wave synthetic diagnostic is lower-36 kHz. The GAM frequency given directly by the code (42 kHz) is also close to the experimental value.

We have also compared the DR signal intensity ´Ps ( f ) df dependence on the vertical antenna displacement, which is expressed in figure 8 in terms of the poloidal wavenumber of the turbulence. As shown in figure 8, the decrease in the DR signal power with growing fluctuation poloidal wavenumber is substantially faster than in the experiment for both synthetic diagnostics. The difference could be attributed to at least two possible reasons. The first reason is related to incorrect modeling of the ETG mode-scale fluctuations leading to their suppression, and absence of the inverse energy cascade in the wavenumber space, which can also result in an underestimation of the turbulence poloidal wavenumber spectrum in the intermediate-scale domain κ ϑ 10 cm -1

. It can occur in spite of sufficiently high density of the GK code grid (poloidal cell size 0.07 cm). This could lead to the faster decay of the synthetic signal power at large antenna displacement. The second reason is related to a possible underestimation by the GK code of the total density fluctuation level. As we show in the next section, this could reduce the role of nonlinear effects: in particular, it could slow down the saturation of the computed DR signal growth with the turbulence level increase at small antenna displacements. The same consequences could result from neglecting in full wave synthetic diagnostic multiple wave reflections between the cutoff and the antenna taking place at the equatorial probing.

RCDR CCF comparison

Inspired by a reasonable agreement of the experimental and synth etic DR spectral data we have performed a comparison of the computed and measured RCDR CCFs using the X-mode DR antenna displaced by +2.0 cm from the equatorial plane at the high-field side of the FT-2 tokamak. The corresponding CCFs (shown in figure 9 as functions of the frequency difference in the master and slave channels (a) and of the separation of turning points (b), which was determined using ray tracing for different probing frequencies, as the difference of minimal radii on the two trajectories) appear to be different, whereas the density fluctuations two-point CCF for frequency f sf i = 300 kHz computed in the measurement region based on the GK data and shown in figure 9(b) is close to the measured RCDR CCF, as was mentioned in [21]. However, this coincidence is just accidental. The CCF provided by the fast synthetic DR is much wider than the experimental one, whereas the full-wave synthetic CCF is situated in between, as shown in figure 9(a). The slow decay of the CCF provided by fast synthetic diagnostics shown in figures 9(a) and (b) is most likely determined by the small-angle-scattering contribution to the DR signal [13], which is not suppressed at the angle of incidence corresponding to y a = +2.0 cm. There are two possible reasons for the drastic difference between the experimental and the fast synthetic RCDR CCFs. The first one is related to the influence of nonlinear effects; in particular, the multiple small-angle scattering, or strong probing wave phase modulation, according to theory [15,16,30], coming into play with the growing turbulence level. The importance of this effect for the DR diagnostics' performance was shown recently at the ASDEX Upgrade [9]. The second reason is related to the GK code spatial grid limitations leading to an underestimation of the turbulence radial wavenumber spectrum width, and thus, to the overestimation of the RCDR CCF width.

The role of the probing wave phase modulation in the FT-2 X-mode RCDR experiment was studied in [31], where the random phase perturbation of the backscattered wave was computed using the geometrical optics approximation and measured. The standard deviation of the computed DR signal phase was shown to be 0.93, whereas for the measured phase it is equal to 1.23. The latter value is close to the computed one and to the level π/3, which indicates the borderline for the DR diagnostics' transition into the nonlinear regime of operation. Strictly speaking, in this parameter domain the reciprocity theorem is not applicable. The obtained agreement of the measured and computed DR signal random phase standard deviation and its relatively high value demonstrates the importance of the first, based upon nonlinear effects, explanation of the drastic disagreement of the experimental and fast synthetic RCDR diagnostics results, and appeals for the application of the full-wave synthetic diagnostics.

As shown in figure 9(a), the CCF provided by the fullwave synthetic diagnostics is much closer to the experimental one. This agreement shows that the second reason mentioned above probably plays a minor role, thus providing an argument in favor of a correct description of the turbulence radial wavenumber spectrum by the GK code with the utilized grid. It should be mentioned that an artificially imposed increase in the density fluctuation level provided by the GK code by a factor of 1.5, roughly corresponding to the excess of the experimental DR signal phase root mean square over the computed one, results in an even better agreement of the experimental and synthetic CCF shown in figure 10. At such a high level of fluctuations the DR already operates in the saturation regime for antenna positions close to the equatorial plane, which leads to the DR signal power growth with the fluctuation level weaker than the quadratic. The same consequences could result from the aforementioned multiple wave reflections neglected in both synthetic diagnostics. The dependence of the DR signal power on the antenna vertical shift in this case is shown in figure 11. Due to the saturation of the DR signal growth close to the equatorial plane the synthetic dependence provided by the full-wave diagnostics in this case appears to be much closer to the experimental one, thus confirming a qualitative explanation for the disagreement given in the previous section. It should also be noted that the increase in the turbulence level by a factor of 1.5 did not alter the DR signal frequency spectra, which remained similar to the exper imental spectra. Thus, the threshold of the strong influence of the strong probing wave phase modulation on the width of the Doppler frequency spectra in the FT-2 case substanti ally exceeds the corresponding thresholds for the RCDR and poloidal wavenumber spectra measurements. This result appears to be in agreement with theoretical expectations [15,30], because the latter thresholds are determined just by transition of reflectometry to the saturation regime, whereas for the former the spectral broadening of the probing wave due to strong phase modulation obtained in [30] should exceed broadening due to other mechanisms. Namely, in the FT-2 experiment it should exceed the broadening due to the strong GAM excitation level, which is dominant.

Conclusions

Based on the results of the paper one should conclude that the ELMFIRE code benchmarking against the FT-2 X-mode DR experimental data, utilizing the fast linear and full-wave DR synthetic diagnostics, has demonstrated good agreement between the measured and computed DR frequency spectra. For all used incidence angles both the spectra frequency shift and width, and in many cases the spectra shapes, were similar. A meaningful difference was only observed at the largest antenna displacement. This difference is probably due to a contribution of small-scale drift modes excited in the plasma at a low level. These modes rotating in the direction of the ion diamagnetic drift are probably not described well by the GK code. Another possible explanation for this difference could be associated with a non-Gaussian form of the antenna beam at the -10 dB level.

ORCID iDs

The mean fluctuation velocities in the DR measurements and in the synthetic diagnostics were close. The GAM characteristics provided by the measurements and by the synthetic DR also appear to be close.

However, the variation in the DR signal power with growing incidence angles in the experiment is slower than that predicted by both of the synthetic diagnostics. This difference could be presumably attributed to the underestimation of the turbulence level and related DR signal saturation effect, and to the overestimation of the turbulence poloidal wavenumber spectrum decay rate with growing wavenumber due to an incorrect description of the ETG mode-scale domain in the GK modeling. The same effect could result from multiple wave reflections between the cutoff and the antenna taking place at the equatorial probing being neglected by the full-wave synthetic diagnostics, as well as by the fast one.

The substantial difference was also found in the decay of the DR signal CCFs with growing frequency shift in the probing wave channels. The quick decrease in the RCDR coherence observed in the experiment and full-wave synthetic diagnostic compared to the fast synthetic RCDR is attributed to nonlinear effects, which manifest themselves in the phase modulation of the probing wave in the former two cases. It was shown to be close to π/3, both in the specially performed measurements and in the GK modeling [31]. In spite of the fact that this value indicates only the beginning of the transition to the fluctuation reflectometry nonlinear regime (and therefore has no impact on the DR spectra) it already has a strong influence on the RCDR performance.

Figure 1 .

 1 Figure 1. The discharge parameter profiles. Density-solid curve; safety factor-dashed curve. Electron temperature-circles; ion temperature-triangles.

Figure 2 .

 2 Figure 2. The DR experiment geometry and the probing wave poloidal electric field distribution computed at y a = 1.5 cm.

Figure 3 .

 3 Figure 3. A comparison of DR spectra for different antenna vertical displacement. Circles-experiment; triangles-fast synthetic DR, squares-full-wave synthetic DR.

1 )Figure 4 .Figure 5 .Figure 6 .

 1456 Figure 4. Dependence of the DR signal frequency shift on the antenna vertical displacement y a and fluctuation wavenumber. Circles-experiment; triangles-fast synthetic DR; squares fullwave synthetic DR.

Figure 7 .Figure 8 .

 78 Figure 7. Power frequency spectra of the velocity oscillations shown in figure 6. Thick line-computed by the GK code, circlesmeasured; triangles-provided by fast synthetic diagnostics; squares-provided by full-wave synthetic diagnostics.

Figure 9 .Figure 10 .

 910 Figure 9. (a) The RCDR CCF against the channel frequency separation. Triangles-fast synthetic DR; squares-full-wave synthetic DR; circles-experiment. (b) The RCDR CCF against the turning point radial separation. Triangles-fast synthetic DR; circles-experiment; double crosses-turbulence CCF by the GK code.

Figure 11 .

 11 Figure 11. Dependence of the BS power on the fluctuation poloidal wavenumber. Circles-experiment; squares-full-wave synthetic DR; stars-FWS DR with the density fluctuation level multiplied by 1.5.

  O.L. Krutkin https://orcid.org/0000-0002-9124-2659 A.B. Altukhov https://orcid.org/0000-0002-5913-6004 A.D. Gurchenko https://orcid.org/0000-0002-2540-2511 T.P. Kiviniemi https://orcid.org/0000-0001-9753-3388 P. Niskala https://orcid.org/0000-0001-7308-9921 C. Lechte https://orcid.org/0000-0003-4171-888X

Acknowledgments

The financial s upport of the Russian S cience Foundation grant 17-12-01110 is acknowledged. The work has been supported by the Academy of Finland grants 278487, 296853 and 318314. CSC-IT Center for Science and Marconi-Fusion HPC are acknowledged for generous allocation of computational resources for this work. The FT-2 tokamak maintenance was supported by the Ioffe Institute, whereas M.A. Irzak acknowledges the support of the Russian Academy of Science Presidium program.


Axel Laureau
email: axel.laureau@cea.fr

Laurent Buiron

Bruno Fontaine

Vincent Pascal

Fission Matrix Interpolation for the TFM Approach Based on a Local Correlated Sampling Technique for Fast Spectrum Heterogeneous Reactors

Reactor safety analysis requires complex coupled multiphysics calculations modeling the full core evolution during various accidental situations. Such calculation codes rely on physical models adapted to the physics of the simulated reactor concept. In this frame, the innovative Transient Fission Matrix (TFM) neutronic approach, based on a conversion of the Monte Carlo response in Green functions characterizing the local transport in the reactor and initialy developed for other kinds of reactors, has been successfully applied on one dimensional sodium fast reactors using the correlated sampling technique as detailed in this article. This method takes into account the sensitivity of the neutron transport to local perturbations in order to estimate on the fly the matrices variations due to a given perturbation distribution. This model provides an accurate estimation of the reactivity and of the power distribution in the core during the accident, including the influence of the sodium density distribution in the plenum. This article also illustrates the applications of this neutronic model on sodium fast reactors for small to large perturbations of the system, and point kinetic feedback coefficients evaluations.

I. INTRODUCTION

To perform reliable reactor safety analysis, accidental transients such as Unprotected Transient Over Power (UTOP) or Loss Of Flow (ULOF) accidents have to be calculated with coupled multiphysics modelling codes at the full core scale. This requires optimized physical models for neutronic and thermalhydraulics at least. In this frame, the innovative Transient Fission Matrix (TFM) neutronic approach has been developed as an hybrid approach between stochastic and deterministic codes, to allow fast and precise calculations during transient simulations. This approach described in [START_REF] Laureau | Transient Fission Matrix: Kinetic calculation and kinetic parameters β eff and Λ eff calculation[END_REF] condenses the spatial and temporal reactor response obtained by Monte Carlo calculations in fission and time matrices that are intended to be used to solve deterministicaly the kinetic equations of the prompt and delayed neutrons during the transient evolution. The reactor behavior is then estimated using a single Monte Carlo calculation prior to the transient evaluation, and finally the neutron kinetics uses directly these matrices without additional calculation. This TFM approach provides a precise estimation of the neutron flux redistribution and of the reactivity during the transient with a reduced computation time.

One of the main challenges is to model the matrix variations during the transient: due to thermal feedbacks the neutron transport is modified. To take this into account without recalculating the matrices at each time step of the transient, a matrix interpolation model has been developed and validated [START_REF] Laureau | Développement de modèles neutroniques pour le couplage thermohydraulique du MSFR et le calcul de paramètres cinétiques effectifs[END_REF][START_REF] Laureau | Coupled neutronics and thermal-hydraulics transient calculations based on a fission matrix approach: application to the Molten Salt Fast Reactor[END_REF][START_REF] Laureau | Transient Fission Matrix: Kinetic calculation and kinetic parameters β eff and Λ eff calculation[END_REF]. This interpolation is appropriate to model Doppler broadening, coolant density or control rods effects for thermal spectrum PWR reactors as well as for the homogeneous Molten Salt Fast Reactor (MSFR). Regarding heterogeneous fast spectrum reactors, the sensitivity to the crossed cells has to be considered on the matrices calculations, due to the large migration area. This has required the development of a new interpolation model.

The present paper describes this new interpolation model based on the calculation of the matrix sensitivity to local feed-back effects, using a correlated sampling technique. Previous work [START_REF] Matthes | Calculation of reactivity perturbations with the Monte Carlo method[END_REF][START_REF] Bernnat | A MONTE CARLO TECHNIQUE FOR LOCAL PERTURBATIONS IN MULTIPLYING SYS-TEMS[END_REF] has shown the efficiency of the correlated sampling technique association to the fission matrices on simple cases. Other work [START_REF] Majumdar | Multiple reactivity calculation using single correlated sampling Monte Carlo simulation[END_REF] illustrated the capability to calculate local density perturbation effect on fission matrices on homogeneous system with one group isotropic Monte Carlo. The model presented here has been developed to extend the TFM approach to the Sodium Fast Reactor (SFR) specificities. This technique consists in a modification of the neutron weight in order to be representative of a perturbed version of the reactor with a single Monte Carlo calculation. The new TFM interpolation model approach, detailed in section II., is based on the use of these perturbed weights to generate the matrices together with their perturbed version. The application of this interpolation model is presented in section III. together with direct Monte Carlo comparisons. The results obtained on the modelling of fuel temperature (Doppler broadening) and coolant (sodium density) perturbations are then detailed. The validity domain of the interpolation for important perturbations such as coolant boiling is also investigated in this paper, comparing the neutron source redistribution and the reactivity variation of the interpolation with direct Monte Carlo calculations. Finally, the interpolation model is also used in section IV. to estimate the distribution of point kinetic feedback coefficients for sodium and temperature perturbations, and also for a clad density perturbation corresponding to a relocalization during accidental situations. These results are presented, together with a comparison to ERANOS [START_REF] Rimpault | The ERANOS code and data system for fast reactor neutronic analyses[END_REF] calculations and a spectral analysis explaining the results obtained.

II. TFM APPROACH AND CORRELATED SAMPLING

Fission Matrices

The fission matrices are the numerical transposition of Green functions aiming to characterise the neutron propagation in the core. The Green function approach consists in reconstructing the neutron propagation as a sum of local contributions, each contribution being the generation system response (νΣ f ψ) to a Dirac-"source neutron emission". Due to the Monte Carlo process, the estimation of these Green functions is based on an arbitrary volume discretization subscripted on i, j, k. For each possible neutron emitted in j, its fission neutron production νΣ f ψ is estimated in each volume i. This production corresponds to the line i -column j term of the fission matrix G that condenses the whole transport neutron information during one generation. The objective of the perturbative approach discussed here is to include the sensitivity to a perturbation localized in volume k on this production.

Note that different kinds of matrices are estimated for the TFM approach. The neutron behavior is different depending on its emission spectrum (prompt or delayed), and the neutron prompt and delayed productions have to be considered distinctly. For this reason, the four matrices

G χ p ν p , G χ p ν d , G χ d ν p ,
G χ d ν d are estimated during the same Monte Carlo calculation, tagging the neutron type for the neutron spectrum and scoring both neutron multiplicities ν p Σ f ψ and ν d Σ f ψ. Note that the neutron spectrum refers to the emission position j and the neutron multiplicity refers to the neutron production position i. Finally the neutron propagation time is also stored in the time propagation matrix. These five matrices are used by the TFM approach to perform neutron kinetics, the fission matrices providing the generational variation of the prompt and delayed neutron distributions, and the time matrix providing the temporal information associated to the generations.

Correlated Sampling Technique

The estimation of a small difference between two different Monte Carlo calculations is suffered from statistical error since the statistical errors are independent. For example for the reactivity, a small perturbation on the sodium density that has an impact of one pcm would require two calculations with a statistical error much smaller than one pcm which is not realistic.

The correlated sampling technique is being used here to improve the estimation of the impact of a perturbation on the neutron path. As mentioned, it consists in associating a perturbed weight to each neutron in order to be representative of a perturbed version of the reactor. This weight is modified at each interaction to take into account that each event doesn't have the same probability to occur in the reference reactor and in a perturbed one. In this way, the different scores can be estimated using a single Monte Carlo calculation for both the reference and the perturbed system with the same neutron histories. Finally, these two estimations have the same statistical error, thus their difference provides a much more precise evaluation of the perturbation on the reactivity or any other scored quantity. Note that if the considered material is not locally perturbed or if the region is voided, the neutron perturbed weight is not modified.

Two kinds of weight perturbation are considered. Considering a given reaction r on the nucleus n at the distance d, with a cross section Σ n,r and a total cross section Σ tot , the probability to sample this distance is Σ tot exp(-d • Σ tot ) and to sample this reaction is Σ n,r/Σ tot . The neutron weight modification due to a perturbation is the ratio of the perturbed probability and the reference one. Considering the distance sampling, the neutron perturbed weight is multiplied by:

Σ pert tot exp(-d • Σ pert tot ) Σ tot exp(-d • Σ tot ) (1)
Considering the reaction type sampling, the neutron perturbed weight is multiplied by:

Σ pert n,r • Σ tot Σ pert tot • Σ n,r (2)
This is a multiplicative process, the neutron weight being modified after each interaction. Each fission has a perturbed weight that results from the multiplication of all the weight modifications of the neutron during its propagation. This perturbed weight is transmitted to the new neutrons produced by the fission, and then to all the future neutron generations, in order to take into account the perturbation of the neutron source. To avoid the accumulation over generations of a large difference between the perturbed weights of the simulated neutrons, this transmission of perturbed weights is limited to a maximum number of generations depending on the reactor.

Correlated Sampling Application to the Fission Matrices

The fission matrices require the estimation of all the i,j elements that correspond to a νΣ f ψ score discretized on the fission position i and neutron creation j volumes. This νΣ f ψ quantity is multiplied by the perturbed neutron weight to generate the matrices of the perturbed reactor.

In order to model the effect of various local perturbations on the fission matrices with the same Monte Carlo calculation, different perturbed weights are associated to each neutron. The perturbed weights are discretized according to each volume k, each associated weight corresponding to a given local perturbation in k. Finally, with two different perturbations (e.g. coolant density and media temperature), the number of perturbed weights associated to each neutron is 2 • k.

All the fission and time matrices are generated for each perturbed weight in k. For the matrix G χ p ν p , we write G den k χ p ν p and G dop k χ p ν p the "variation matrices" generated with the correlated sampling technique, corresponding to the variation of G χ p ν p due to a local perturbation of the coolant density and of the temperature in k.

The neutron weight propagation through generations, required for the correlated sampling technique, has been implemented and can be used for the calculation of G den k χ p ν p and G dop k χ p ν p . However, thanks to the capability of the fission matrices to predict the perturbed neutron source through the Eigen vector, this propagation of the neutron weight is not required if the mesh discretization is fine enough and the generated matrices are identical with and without neutron weight propagation. For a given configuration, a sensitivity study on the mesh discretization is required.

Interpolation of the Fission Matrices

In order to estimate the modification of the fission matrices due to an arbitrary perturbation amplitude and distribution, the variation matrices can be used to interpolate the fission matrices. The global perturbation is reconstructed as a sum of local perturbations, together with a linear interpolation for the density effect and a logarithmic interpolation for the Doppler effect. For any matrix G χ x ν x of the TFM approach, the interpolated matrix with a density and a temperature perturbations is calculated using Eq. [START_REF] Majumdar | Multiple reactivity calculation using single correlated sampling Monte Carlo simulation[END_REF].

G χ x ν x (∆ρ sodium (k), T(k)) = G χ x ν x + k G den k χ x ν x • ∆ρ sodium (k) + k G dop k χ x ν x log( T (k) /T ref (k)) log( (T ref (k)+300) /T ref (k)) (3)

Source Neutron Redistribution and Reactivity Estimation

The interpolated matrix is representative of the local neutron propagation in the perturbed system. Considering a time step t associated to a perturbation shape (∆ρ sodium (k), T(k)), the Eigen vector N p (t) of the matrix

G χ p ν p (t) = G χ p ν p (∆ρ sodium (k), T(k)
) is solution of the Eq. ( 4) and corresponds to the equilibrium neutron distribution. The associated Eigen value is the prompt multiplication factor of the system.

G χ p ν p (t)N p (t) = k p (t)N p (t) (4)
In order to take into account the coupled effects between prompt and delayed neutrons, we use the G all (t) block-matrix (Eq. 5) constructed with the "prompt to prompt" and "delayed to delayed" matrices on the diagonal, and the crossed productions out of the diagonal. The Eigen vector N eq,all (t) = N p,eq (t) N d,eq (t) corresponds to the prompt and delayed distribution in the reactor at equilibrium, and the associated Eigen value is the effective multiplication factor k eff (t) of the reactor.

G all (t) =          G χ p ν p (t) G χ d ν p (t) G χ p ν d (t) G χ d ν d (t)          (5)

III. APPLICATION OF THE TFM APPROACH WITH THE PERTURBATIVE INTERPOLATION MODEL 1. Application Case Presentation

The development of this perturbative interpolation model has been performed in the general framework of the low void effect [START_REF] Sciora | Low void effect core design applied on 2400 MWth SFR reactor[END_REF] sodium fast reactor ASTRID. This kind of concept is characterized by a large sodium plenum increasing the neutron leakage in case of a sodium density decrease.

A simplified geometry of a 1D assembly has been considered as application case to check the capability of this approach Fiss -1500 K Na -600 K B 4 C -600 K 16 to be representative of the physics involved in this kind of reactor with the variation of the flux gradient in case of sodium loss. The geometry, composed of three areas (fissile, sodium and B 4 C), is represented in Fig. 1 and the material temperatures and isotopic compositions are given in Table I. These compositions are considered radially homogeneous. The geometry boundary condition is a radial reflexion and axial leakage. Note that the fuel region here corresponds to an assembly homogenization so that it contains fuel, sodium and steel. All the calculations presented in this paper have been performed with a modified version of the neutron Monte Carlo code Serpent [START_REF] Leppänen | The Serpent Monte Carlo code: Status, development and applications in 2013[END_REF] and with the JEFF3.1 database. One billion of neutrons have been simulated per calculation and a 60 bins geometry discretization for i j and k is considered. A sensitivity study has been performed: the results detailed in the article are identical using a finer 120 bins mesh, and also using a propagation of the neutron perturbed weight on 8 generations. The source neutron distribution in the system is presented in Fig. 2 (left). We can notice the asymmetry of this distribution in the fuel due to the large neutron leakage at the bottom of the geometry (left part of the figure). The neutron reflection contribution of the sodium plenum is directly visible with the large amount of neutrons created in the fuel near to the sodium: without the reflection of the sodium, the neutron source at 50 cm would be equivalent to the one at 0 cm. For this reason, the leakages are increased accordingly if the sodium density is reduced.

The neutron production shape corresponds to the Eigen vector of the fission matrix G χ p ν p displayed in Fig. 2

(right).

Each column j of the fission matrix corresponds to the propagation of the neutron in the i volumes. Considering an emission in the center of the fuel (ie. column j = 15 of the matrix), each line of this column represents the probability that a neutron produces a fission neutron during one generation. This proba-bility is zero in the sodium (bins 30-44) and in the B 4 C areas (bins 45-59). The values on the diagonal (j = i) correspond to the neutron production at the same place as the neutron emission, which is the most probable fission position. The general diagonal shape of this matrix shows the propagation of the fissions around the neutron emission, the values far from the diagonal correspond to large migration areas.

Matrix Calculations

Using the correlated sampling technique, different variation matrices are generated. The matrices presented in Fig. 3 are the variation of the G χ p ν p matrix due to a local perturbation in k of the sodium density or of the fuel temperature. We consider in this paper a reference estimation of these matrices for a sodium density perturbation of -1% and for a temperature perturbation of +300K. A local variation of the sodium density of -1% produces a redistribution of the neutron propagation with a strong discontinuity around the perturbation position k. If the emission and the fission position are on the same side, (i and j) > k or (i and j)< k, the impact of the perturbation is a decrease of the fission neutron production. At the opposite, if the perturbation is located between i and j, the neutron production is increased because the spectrum becomes harder. The effect of a sodium density decrease in the sodium plenum is a strong fission production reduction near to the sodium area due to the neutron leakage in the B 4 C.

Concerning the Doppler effect, it produces a global decrease of the fission production around position k. The Doppler broadening favours the absorptions, and this effect is more important for a perturbation near to the sodium area where the spectrum is softer.

Interpolation of the Matrices and Validation

In order to validate the interpolation model and the implementation of the approach, two configurations have been considered. A perturbation of the reactor with a different shape and amplitude has been calculated with a variation of the sodium density of -2% in the sodium area, and a variation of -300K at the top of fuel area. In this way this checks the interpolation based on the reference variation of the sodium density of -1% and of the temperature of +300K (reference matrices). The interpolated fission matrices corresponding to these perturbations are calculated using:

G χ x ν x -2%Na = G χ x ν x + k∈{30,44} G den k χ x ν x • 2 G χ x ν x -300K = G χ x ν x + k∈{15,29} G dop k χ x ν x log( 1200 /1500) log( 1800 /1500) (6)
We remind here that the reference matrices are calculated using direct Monte Carlo calculations, without using the correlated sampling technique. Two distinct Monte Carlo calculations are performed, for the perturbed and the reference system, and the difference of the generated fission matrices provides the reference values for G χ x ν x -2%Na and G χ x ν x -300K .

The results obtained are presented in Fig. 4. The results of the interpolation model are represented on the left, the reference estimation using direct Monte Carlo calculation are represented in the middle, and the difference between these two estimations are represented on the right. For both density and Doppler effects, the global shapes are very well reproduced by the interpolated matrices and the differences are due to the statistical errors only.

In addition to the fission matrix comparison, a second point is the validation of the redistribution shape prediction. For this, the Eigen vector of the interpolated matrix is compared to the estimation of the source neutron shape difference of the two Monte Carlo direct calculations. The two configurations are represented in Fig. 5. For the density perturbation (red) and the Doppler perturbation (blue), the interpolated results (dashed line) perfectly reproduce the Monte Carlo calculation (solid line).

Finally, the Eigen values associated to the interpolated matrices are compared to the reactivity variations of the Monte Carlo direct estimation in table II. Once again the interpolation model shows its capability to precisely predict the reactivity variation due to a perturbation in the system.

Interpolation of the Matrices for Large Perturbations

The fission matrice interpolation uses variation matrices G χ x ν x calculated with a small perturbation around the nominal state. Using these matrices may induce a bias if we consider large perturbations of the reactor compared to the nominal state due to a potential non-linear behavior. Note that this interpolation model does not assume that the flux redistribution is linear at the system scale, but that the local system response on one generation is linear. In this section, large sodium variations are considered as illustration since this phenomenon produces the larger flux in the core.

Two characteristic amplitudes are considered: a variation of -15 % that is the order of magnitude of the sodium density variation before sodium boiling, and a variation of -100 % that corresponds to a total void effect. Four perturbation lo-Case density 50-75 cm Doppler 25-50 cm ∆ρ ref -60 ± 3 69 ± 3 ∆ρ interpolation -60.7 67.8 difference (1 ± 5) % (-2 ± 5) % TABLE II. Reactivity variation due to a local variation of the sodium density and of the temperature calizations have been calculated in order to obtain distinct characteristic behaviors: case A in the bottom of the fissile region (0 cm up to 25 cm), case B in the upper part of the fissile region (25 up to 50 cm), case C in the first 10 cm of the sodium region (50 cm up to 60 cm), and case BC that combines B and C (25 cm up to 60 cm). Each case represented in Fig. 6 is declined for the two amplitudes mentionned. The results obtained with the sodium density perturbation at -15 % and -100 % are respectively presented in Figs 7 and8, and the reactivity in Tables III and IV. As previously illustrated, the non-zero terms of the matrices are only in the 0-30 bins (corresponding to the fissile area in 0-50 cm). For each case, the reference calculation (middle) corresponds to two direct Monte Carlo calculation with the corresponding sodium modification shape and amplitude: the variation matrix is the difference between the perturbed and the reference matrices. Note the two distinct Monte Carlo used for this reference calculation are performed with one billion of neutron, while all the matrices required for the interpolation (left) are calculated using only one calculation with also one billion of neutrons.

Different observations can be done concerning the physics of the matrices variation independently of the perturbation amplitude (Figs. 7 and8):

• Concerning case B, a decrease of the fission probability along the diagonal is observed in the upper fissile area. This effect is due to the increase of the neutron mean free path. An important increase of the fission probability is observed far from the diagonal due to the sodium density decrease that implies a harder neutron spectrum.

• If the sodium perturbation is near to the bottom of the fissile area (case A), the effect on the matrix is almost the same but the proximity of the bottom boundary reduces the fission probability far from the diagonal.

• For a perturbation in the sodium area (case C), the effect is a depletion of the fission matrix near to the sodium.

• For a perturbation combining cases B and C, the resulting matrix is the sum of the two corresponding matrices. Case III. Reactivity variation due to -15 % sodium variation from TFM and the Monte Carlo reference calculation A perturbation in the sodium area implies an important decrease of the neutron source redistribution near to the upper fissile area (Fig. 7 -right, case C ) due to the large increase of the neutron leakages in the B 4 C. For a perturbation in the fuel, two opposite effects occur: the local mean free path increase tends to locally decrease the neutron redistribution, while a harder neutron spectrum tends to increase the fission rate. For this reason, we can see on cases A and B that the redistribution is almost negative in the perturbation position, and positive at the boundaries of the perturbation with a large influence of the material next to the perturbation (sodium, unperturbed fuel or leakage).

∆ρ ref ∆ρ interpolation difference A -91 ± 2 -95 4 ± 2 % B 92 ± 2 88 -4 ± 2 % C -181 ± 2 -183 1 ± 1 % BC -94 ± 2 -96 2 ± 2 % TABLE
With the -15 % sodium density variation, the amplitude of the variation matrices is limited compared to the -100 %, and thus the reference and the perturbed fission matrices estimated using two distinct Monte Carlo calculations are not significantly different. Since the statistical errors are independent, Fig. 8. TFM interpolated matrix variation (bins 0-30) associated to the -100 % sodium variation (left), reference matrix variation (middle) and comparison of the associated source neutron redistributions. IV. Reactivity variation due to -100 % sodium variation from TFM and the Monte Carlo reference calculation the reference variation matrix that is the difference between the two calculated matrices has a large statistical error. This effect is reduced for the -100 % since the residual statistical error is much smaller than the difference of the perturbed and reference matrices. A good global agreement is obtained in both cases, the behavior of the neutron local propagation is correctly reproduced. Note that the feedback is positive for case B since the harder spectrum tends to increase the reactivity, while it is negative for case A due to the increased neutron leakages. This effect will be discussed in section IV. with the point kinetic feedback coefficients estimation.

Case

The neutron source redistributions (right column) are also correctly reproduced. The agreement is very good for a -15 % perturbation, and a small difference appears for the total void configuration.

The reactivity variations are also reproduced with a very good agreement for the -15 % perturbation (Table III), the maximum difference being of 4% on the reactivity variation. The limit of the interpolation appears with a -100 % perturbation (Table IV). For a total void effect, the reactivity variation is predicted with a maximum bias of 16 %. Even if the global shape of the matrices is good and then their Eigen vectors (source neutron distributions) are correct, the Eigen values (multiplication factors) are very sensitive to the non-linearity for a local total void configuration. A possible solution that can be investigated for future developments is to estimate the variation matrices at different intermediate densities between 100% and 0%, and finally to interpol the final matrix with these intermediate matrices.

IV. POINT KINETIC LOCAL FEEDBACK DISTRIBU-TION 1. Feedback Calculations

The objective of the initial development of the TFM approach was the calculation of effective kinetic parameters and the capability to perform spatial kinetic calculations. From the new developments presented here, local feedback coefficients can also be calculated by using the variation matrices. The TFM approach can thus provide Monte Carlo based point kinetics precise parameters to other calculation codes.

Each variation matrix corresponds to a local perturbation in the reactor (material density and Doppler broadening in this work). The effective multiplication factor associated to this perturbation is directly given by the Eigen value of G all + G k all . Finally the reactivity variation associated to each local perturbation in the volumes k is calculated using TFM.

Note that the corresponding Eigen vector can also be calculated for each k value. Assuming that the impact of the variation of the neutron source distribution due to a system modification in cell k is negligible on the propagation of another modification in cell k in the reactor, then a pseudo spatial kinetic model can be established. This model uses at each time step a sum of variation vectors weighted by the amplitude of the system modification in k. In this way the calculation time would be very fast, but this approach can not be able to observe neutron decoupling effect.

In order to have a comparison for the feedback coefficients obtained, S N calculations with 33 energy groups have been performed with the ERANOS code [START_REF] Rimpault | The ERANOS code and data system for fast reactor neutronic analyses[END_REF]. These deterministic calculations are based on the traditional two level lattice/core scheme and the JEFF 3.1 nuclear database. First, self-shielded cross sections are computed by the ECCO code cell, using the fundamental mode assumption for all the materials of the 1D core geometry. For fissile materials, a buckling search algorithm is used to obtain the critical flux for the cross section collapsing to a 33 energy group mesh. For the subcritical materials such as fertile or structural parts of the 1D sub-assembly, the process is based on source calculations using the spectrum coming from previous fissile calculations. The boundary conditions are a flux leakage on the axial boundaries for both cases, and no volumetric leakage are considered since there is no radial leakage for this application case.

Figure 9 shows the density (left) and Doppler (right) feedbacks coefficient calculated for the system. The density feedbacks are calculated for a sodium and clad density variation of -1 %. This Doppler feedback is calculated with a logarithmic dependency due to the temperature perturbation amplitude in the correlated sampling process (+300 K). A very good agreement is obtained between TFM (red) and ERANOS (blue). Concerning the density effect, the feedbacks are different in the three regions. In the fuel region the sodium (plain line) and the steel (dashed line) feedbacks have the same behavior, the amplitude being larger for the steel due to the larger absorptions in this media. Near to the fuel boundaries, the feedback is negative since a decrease of the nuclei density increases the local neutron mean free path, and then the leakages. In the middle of the fuel region the feedbacks are positive, indeed the decrease of the sodium density makes the neutron spectrum harder, together with the absorption decrease in the steel if its density is reduced.

In the B 4 C, the feedbacks are negative since the presence of sodium and steel reduces the absorptions on the boron due to the reflections. In the sodium plenum, the only component is sodium without any steel and the clad feedback is thus null here. We can see that the sodium feedback value is constant in this region. The impact on the reactivity is exactly the same if the perturbation is near to the fissile (eg. at 55 cm) or near to the B 4 C (eg. 70 cm) even if the neutron flux and importance are reduced in the latter. Such an effect was expected with this one dimensional geometry without radial leakages. The only leakages are axial, and for the reason, the probability that a neutron leaving the fissile area arrives in the B 4 C only depends on the total amount of sodium nuclei between these two regions. The reflector capability of the sodium in this region only depends on the total amount of sodium and not of its distribution.

Concerning the Doppler effect, this effect is mainly negative since the Doppler broadening favors the absorption/fission ratio. Between 0 cm and 5 cm the Doppler effect is positive: the total cross section increase limits the leakage on this part of the geometry. An important local variation may be observed near to the fissile-sodium interface. Indeed the reflected neutrons returning from the sodium area have a softer spectrum which directly impacts the Doppler broadening sensitive to the neutron spectrum.

A slight difference is obtained between TFM (red) and ER-ANOS (blue). This difference is due to the important spectral variation detailed in the next section. The neutron spectrum variation is also large between the sodium and the fissile areas, meaning that the fundamental mode assumption may no be appropriate. Finally, another possible origin for this discrepancy is the self-shielded multigroup cross section anisotropy. A neutron reflected from the sodium down to the fuel area comes back with a smoother spectrum, and due to the progressive spectrum variation inside the fissile, the self-shielded cross sections may be different depending on the neutron direction.

Spectral analysis

In order to confirm the origin of the difference between deterministic and TFM based Doppler feedback coefficients, spectrum maps have been generated and are detailed here. These maps, evaluated with the Serpent Monte Carlo code, consist in a flux score in the system with a double discretization: axially (500 bins in ordinate) and in lethargy (5000 bins in abscissa). The Serpent code has been modified to obtain an additional information for each neutron: the direction of its axial velocity, positive if the neutron goes from the bottom to Fig. 10. Neutron flux normalized per source neutron: the "A" map represents the total flux, the "B" and "C" maps represent the neutrons with respectively a positive (up) and a negative (down) axial velocity, and the difference of "B-C" is presented in map "D" in percents of the average values. the top of the assembly, and negative elsewhere.

The first maps (Fig. 10) represent the flux repartition in the system. The A map corresponds to the direct flux for each lethargy-position bin, while the B and C maps represent the map flux of the neutrons with respectively a positive and negative axial velocity. The D map represents the relative difference between B and C.

We can see that the mean neutron energy in the fissile region (between 0 and 50 cm) is equal to around 300 keV even if the neutrons are produced here at around 2 MeV, due to scattering in the fuel. Then the neutrons that move to the sodium area are thermalized as observed between 50 and 70 cm. The fast component (1-10 MeV) is cut within a few centimeters in the sodium. The component around 100 keV is more important and we can see the accumulation of a neutron population around 1keV, below the sodium absorption resonance that depletes the neutron spectrum map. We can note the large number of resonances in the fissile area compared to the sodium area, especially at low energy.

As explained, map D represents the difference between the up and down velocity components. Using the S N transport solver for the ERANOS code, the strong flux anisotropy that is directly visible on this D map is correctly taken into account.

We need to go further to understand the difference between ERANOS and TFM. Figure 11 represents the spectrum maps: it corresponds to the flux map but normalized in each position bin in order to see how the neutron spectra are modified inside the geometry. Fig. 11. Neutron spectra normalized per axial bin: the "A" map represents the total spectra, the "B" and "C" maps represent the neutrons with respectively a positive (up) and a negative (down) axial velocity, and the difference of "B-C" is presented in map "D" in percents of the average values.

We can see on map D that the neutron spectrum needs a few tens centimeters to be converged on a "fundamental spectrum" between each material interface. Additionally, we can see a strong spectral anisotropy in the energy range of the resonances. For these reasons, the fundamental mode assumption may be not appropriate on this kind of geometry where the fuel height is small, and new calculation schemes considering anisotropy and axial dependence in the self-shielded cross sections should provide better results.

V. CONCLUSIONS

The present work provides an innovative way to predict the neutron redistribution and the reactivity variation due to a perturbation in a fast spectrum heterogeneous reactor. The paper demonstrates the capability to interpolate fission matrices on the fly for a given shape and amplitude of coolant density and material temperature using a sum of local contribution. Thanks to a correlated sampling technique, the effect of different local perturbations on the fission matrices are estimated using a single Monte Carlo calculation and with a small statistical error. The precision obtained with this model on the neutron redistribution and on the reactivity is around the percent on a partial assembly representative of the low void effect of the ASTRID reactor.

Perturbations with large amplitude have also been tested to investigate the limits of this interpolation model. A sodium variation up to -15 % that corresponds to the limit before boiling is reproduced with a very good agreement compared to direct Monte Carlo calculations. A variation of -100 % that corresponds to a local total void effect is predicted with a good agreement on the neutron redistribution but a bias of 16 % on the reactivity variation due to the non-linearities.

Finally point kinetics feedback distributions have been sucessfully calculated with this Monte Carlo based approach, and compared to the ERANOS deterministic code.

This model can be used for future applications to perform precise neutron kinetics with a coupling to the thermalhydraulics, only one Monte Carlo calculation being required prior to the transient and all the interpolations being performed with a reduced computation time.

Fig. 1 .

 1 Fig. 1. Case geometry description in cm.

Fig. 2 .

 2 Fig. 2. Source neutron distribution (left) and reference G χ p ν p matrix (right).

Fig. 3 .

 3 Fig. 3. Variation matrices G den k χ x ν x (top) and G dop k χ x ν x (bottom) for a perturbation in volumes 15 (left), 25 (middle) and 40 (right).

Fig. 4 .

 4 Fig. 4. Fission matrix variation due to a local variation of the sodium density of -2% in the sodium plenum (top) and to a local variation of the temperature of -300K (bottom): results of the interpolation model (left), direct Monte Carlo calculation (middle) and difference (right).

Fig. 5 .

 5 Fig. 5. Source neutron redistribution due to a local variation of the sodium density (red) and of the temperature (blue): results of the interpolation model (dashed line) and of the direct Monte Carlo calculation together with the statistical error (solid line).

Fig. 6 .

 6 Fig. 6. Position of the perturbation in cases A,B,C and BC.

Fig. 7 .

 7 Fig. 7. TFM interpolated matrix variation (bins 0-30) associated to the -15 % sodium variation (left), reference matrix variation (middle) and comparison of the associated source neutron redistributions.

Fig. 9 .

 9 Fig. 9. Density feedback coefficient (left) associated to a local perturbation of -1 % in the sodium (plain line) and the steel (dashed line), and Doppler logarithmic coefficient (right) calculated as ∆ρ log( T +300 /T) , with TFM (red) and ERANOS (blue).

VI. ACKNOWLEDGMENTS

The authors wish to thank the CEA (Commissariat à l'Energie Atomique et aux Energies Alternatives) and the IN2P3 department of the CNRS (National Center of for Scientific Research) for their support.


C Le Gall

M Cologna

M Holzhäuser

K Popa

J Boshoven

H Hein

J.-L Hazemann

O Proux

C Riglet Martial

J Léchelle

F Audubert

Y Pontillon

An innovative Method to study Volatile Fission Products speciation in nuclear fuels under Severe Accident Conditions

Keywords: Severe accident, Fission products, volatile, Cesium, SPS, in-situ XAS

The development of realistic models for fission products behavior during a severe accident requires experimental data on fission products speciation into the fuel (as irradiated and during the transient) together with thermodynamic predictions. Due to the limitations in terms of experiments and characterization techniques available to study fission products speciation in nuclear irradiated fuels, simulant materials (SIMFUEL manufactured thanks to conventional process) are often used. However, these materials do not allow the study of volatile fission products because they are totally released during the sintering stage. In the case of SIMFUEL with implanted fission products surrogates, the chemical state and the corresponding evolution of the phases during the thermal sequence are not entirely representative of the real case. An innovative method to study volatile fission products speciation in SIMFUEL samples has thus been developed. Spark plasma sintering was used to synthesize dense SIMFUEL samples containing cesium in collaboration with the JRC-Karlsruhe within the frame of the GENTLE program. Thermodynamic calculations made using the TAF-ID v.6 associated with the FACTSage software enabled to determine the conditions in which to perform thermal treatments. Intemperature XAS experiments in conditions representative of a sevre accident are developed in collaboration with FAME-UHD beamline at the ESRF (Grenoble, France). Additional thermal treatments will also be performed to enable further characterizations (ceramography, SEM-EDX…).

Introduction

During irradiation cycles of nuclear fuel in a reactor, Fission Products (FP) are created and retained in the fuel, modifying its microstructure and properties (thermal conductivity, creep behavior…). In case the three containment barriers would fail, as it might be the case during a Severe Accident (SA), FP can be released from the fuel into the environment. Their chemical nature and the quantities released (source term) at these stages strongly depend on the scenario and the accidental conditions involved (temperature, type of fuel, burn-up, oxygen potential…). A deep understanding of FP release mechanisms out of the fuel is thus essential to get a realistic estimation of the source term and to be able to better predict the consequences of the accident.

Therefore, many research programs on irradiated fuels have started worldwide, leading to a better understanding of FP behavior during temperature transients under a controlled atmosphere. A classification of FP according to their volatility has also been proposed largely thanks to the HEVA/VERCORS programs, launched by the CEA in 1989 [START_REF] Pontillon | Behaviour of fission products under severe PWR accident conditions: The VERCORS experimental programme-Part 2: Release and transport of fission gases and volatile fission products[END_REF]:

• Volatile (Xe, Kr, I, Cs, Sb, Te, Cd, Rb, Ag); • Semi-volatile (Mo, Rh, Ba, Pd, Tc); • Low-volatile (Ru, Ce, Sr, Y, Eu, Nb, La); • Non-volatile (Zr, Nd, Pr).

To go further, getting experimental data on FP speciation and confronting them to thermodynamics predictions are necessary. Complementary to irradiated fuels studies, the use of SIMFUEL makes it possible to investigate FP speciation by means of powerful techniques such as X-ray absorption spectroscopy [START_REF] Lucuta | Microstructural features of SIMFUEL -Simulated high-burnup UO2-based nuclear fuel[END_REF]- [START_REF] Geiger | Fission products behaviour in UO 2 submitted to nuclear severe accident conditions[END_REF]. A model based on experimental observations performed on irradiated fuels, SIMFUEL and thermodynamic calculations has thus been proposed to explain the released fractions during a VERCORS-like sequence [START_REF] Geiger | Insights on fission products behaviour in nuclear severe accident conditions by X-ray absorption spectroscopy[END_REF], [START_REF] Kleykamp | The chemical state of the fission products in oxide fuels[END_REF]- [START_REF] Geiger | Study of fission products (Cs, Ba, Mo, Ru) behaviour in irradiated and simulated nuclear fuels during severe accidents using X-ray absorption spectroscopy, SIMS and EPMA[END_REF]. According to this model, reactions between FP and other structural materials take place at different temperatures, strongly influenced by the oxygen partial pressure.

Among the volatile FP, Cs is of great interest for it has a high fission yield and might present a huge radiological impact on the environment and on population on the long term. It is also able to interact with other FP and change their release kinetic. Indeed, it can react with Mo and U to form cesium molybdates or uranates and can also be found in oxide precipitates along with Ba and Zr [START_REF] Kleykamp | The chemical state of the fission products in oxide fuels[END_REF], [START_REF] Kleykamp | Composition and structure of fission product precipitates in irradiated oxide fuels: Correlation with phase studies in the Mo-Ru-Rh-Pd and BaO-UO2-ZrO2-MoO2 Systems[END_REF], [START_REF] Paschol | Phase equilibria in the pseudoquaternary BaO-UO2-ZrO2-MoO2 system[END_REF]. However, some major uncertainties still remain notably concerning Cs behavior. Indeed, the conventional way to synthesize SIMFUEL does not allow volatile FPs to remain in the UO2 matrix as they are released during the sintering process (1650°C, 4h). In the case of implantation processes, the evolution of the phases during the thermal sequence is not representatively reproduced [START_REF] Geiger | Study of fission products (Cs, Ba, Mo, Ru) behaviour in irradiated and simulated nuclear fuels during severe accidents using X-ray absorption spectroscopy, SIMS and EPMA[END_REF], [START_REF] Hocking | Study of fission-product segregation in used CANDU fuel by X-ray photoelectron spectroscopy (XPS) II[END_REF], [START_REF] Hocking | Migration behaviour of iodine in nuclear fuel[END_REF].

Recently, SPS was successfully used to synthesize dense UO2 pellets containing CsI [START_REF] Wangle | Simulated UO2 fuel containing CsI by spark plasma sintering[END_REF]. The objective of the present work is to investigate the possibility of studying Cs behavior inside more representative samples towards irradiated fuels. To this end, an innovative approach which will be described in this paper has been developed and lies on:

• The syntheses of SIMFUEL samples containing Cs, Mo, Ba and Zr through SPS at the Joint European Commission (JRC) Karlsruhe, Germany; • Thermodynamic calculation for dimensioning the thermal treatments;

• In temperature X-ray Absorption Spectroscopy performed under controlled atmosphere on FAME-UHD beamline at the European Synchrotron Radiation Facility (ESRF); • Additional thermal treatments performed into a dedicated furnace in the same conditions in order to verify that post-tests observations are consistent with FP states at high temperatures.

Experimental

Materials

In order to study Cs behavior in SA-like conditions, five batches of samples containing FP surrogates under different chemical forms have been prepared (Table 1). The ratios between each FP are representative of a pressurized water reactor fuel with a burn-up of 76 GWd.tHM -1 . Only the results obtained on batch n°4 will be detailed. Commercial depleted UO2 was first pre-reduced at 800°C during 4h under Ar + 6.5% H2 in order to avoid a stoichiometry gradient to appear in the pellet after sintering (Figure 1) [START_REF] Cologna | On the Role of the Electrical Field in Spark Plasma Sintering of UO2+x[END_REF]. The lattice parameter was determined through XRD on a Rigaken MiniFlux 600 device and the final O/M was calculated according to equation ( 1) [START_REF] Teske | Investigation of the oxygen activity of oxide fuels and fuel-fission product systems by solid electrolyte techniques. Part I: Qualification and limitations of the method[END_REF]. It was reduced from 2.16 to 2.01.

𝑎𝑎 = 5.4705 -0.132𝑥𝑥 (1)
Commercial MoO3 (99.5% purity, MERCK, Figure 2) was used as received and cesium uranate was synthesized according to the protocol described in [START_REF] Osborne | Cesium uranate, Cs2UO4: heat capacity and entropy from 5 to 350 K and standard Gibbs energy of formation at 298.15 K[END_REF]. CsNO3 (MERK) and depleted UO2 (Cogema) were mixed by grinding in an agate mortar and placed in a calcination furnace. The temperature was increased gradually and held at 700°C under air during 24h. The transformation of CsNO3 was not complete at the end of the sequence (Figure 3) so the heating process was repeated a second time up to 800°C. A Rietveld refinement of the XRD patterns showed that a mix between Cs2UO4 (22%) and Cs2U2O7 (75%) was obtained.

Syntheses of SIMFUEL samples

SPS is an electrical field assisted sintering technique that allows obtaining dense UO2 pellets at lower temperature in a shorter time [START_REF] Tyrpekl | Implementation of a spark plasma sintering facility in a hermetic glovebox for compaction of toxic, radiotoxic, and air sensitive materials[END_REF]. This enables the incorporation of volatile FP surrogates in the UO2 matrix and has been deployed at the JRC-Karlsruhe to synthesize SIMFUEL samples [START_REF] Wangle | Simulated UO2 fuel containing CsI by spark plasma sintering[END_REF].

Sample preparation was carried in a glovebox under Ar atmosphere. The additives powders were first ground separately in an agate mortar. The pre-reduced UO2 was then added to the mixture and ground again. The powder was finally poured in a 6 mm SPS matrix containing a graphite foil. A pre-compaction was performed at 500 N (~17.7 MPa) and the following cycle was run:

• Evacuate to <10 Pa • Press to 2.

Characterization of SIMFUEL samples

The first SEM observations were performed on a Tescan Vega LSH device equipped with an Oxford EDX detector. Fractured samples were used in order to get rid of the graphite layer on their surface. Acceleration tension was set to 20 kV. Semi-quantitative analyses have been performed for U, O, Mo and Cs.

Results and discussion

Chemical characterization

Image analysis was performed using the analySIS pro software v.5.0 from Olympus Soft Imaging Solutions. Polyhedral grains with a size from 1 to 4 µm can be seen. The average grain size was measured on 30 grains to be 2.11 ± 0.89 µm (Figure 4). A lot of pores are still visible between the grains. The observed microstructure matches the approximated density calculated for this sample (86.1± 4.3 %d 𝑡𝑡ℎ ) and confirms that the sintering was uncomplete. Source Term issues and the Paper N° (order as in the announcement) Some grains containing Cs, which could be urinate precipitates, have been found in the samples (Figure 4). The solubility of Cs in UO2 measured to be below 0.06wt% [START_REF] Kleykamp | The solubility of selected fission products in UO2 and (U, Pu)O2[END_REF], [START_REF] Bagger | Observations on the Release of Caesium from UO2 Fuel[END_REF] could explain the formation of Cs precipitates.

Figure 4: SEM-SE image of fractured surface of a sample from batch n°4. A grain containing Cs is indicated in red.

Melted area containing Mo can also be observed (Figure 5). This could be explained by the melting of MoO3 occurring at 795°C which is well below the maximum temperature of 1200°C reached during the sintering. The atmosphere inside the SPS furnace might not have enabled the reduction of MoO3 into MoO2 or Mo metal (melting at 1100°C and 2617°C respectively). Furthermore, the solubility of Mo in UO2 is below 0.002 wt% [START_REF] Kleykamp | The solubility of selected fission products in UO2 and (U, Pu)O2[END_REF], [START_REF] Giacchetti | Behavior of Molybdenum in Mixed-Oxide Fuel[END_REF]. The presence of Mo melted phases could thus have hindered the sintering to take place. SEM-EDX analyses and X-ray cartographies on polished samples have to be performed in order to determine accurately the initial FP state. Moreover, chemical analyses by mean of Inductive Coupled Plasma Mass Spectrometry (ICP-MS) will enable to determine the amount of Cs released during the sintering.

Thermodynamic predictions

The conditions in which the thermal treatments will be performed have been determined thanks to predominance diagrams calculated on the FACTSAGE software coupled with the TAF-ID [START_REF]OECD NEA/NSC: Thermodynamics of Advanced Fuels -International Database (TAF-ID)[END_REF]. These diagrams represent the oxygen potential vs the temperature (Figure 6). The brown lines correspond to the limits of the different domains. The ratio between H2O/H2 and CO/CO2 are drawn as the red and blue dashed lines respectively. The black dashed line corresponds to the iso-concentration ratio between Cs2MoO4 and Cs in a logarithmic scale.

The objective of the thermal treatments is thus to probe the different domains of the diagrams to bring direct experimental evidence for the existence of the phases involving the FP of interest.

In-temperature X-ray Absorption Spectroscopy

An in-temperature XAS experiment was developed in collaboration with the beamline FAME-UHD (ESRF, France). The set-up designed to perform the experiment has been extensively described in [START_REF] Richet | Raman spectroscopy, X-ray diffraction and phase relationship determinations with a versatile heating cell for measurements up to 3600 K (or 2700 K in air)[END_REF], [START_REF] Neuville | In situ High-Temperature Experiments[END_REF] and used to carry out in-situ XAS experiments on radioactive materials at the INE beamline (ANKA synchrotron, Germany) [START_REF] Caisso | In situ characterization of uranium and americium oxide solid solution formation for CRMP process: First combination of in situ XRD and XANES measurements[END_REF].

The sample holder has been modified in order to handle and confine a 1 mg sample on a W-5% Re wire heated by Joule effect (Figure 7). A support has been added to avoid deformation of the wire at high temperature. The temperature is monitored during the test by both a thermocouple located beneath the sample and a pyrometer measuring the temperature at the surface of the sample. The atmosphere inside the furnace is controlled thanks to flowing gas mixtures (Ar, H2O/H2 or CO/ CO2). Analyzer Spectrometer (CAS) equipped with five Si(311) crystals will be used as it enables to probe very diluted elements in a heavy matrix [START_REF] Hazemann | High Resolution Spectroscopy on an X-ray Absorption Beamline[END_REF], [START_REF] Llorens | High energy resolution five-crystal spectrometer for high quality fluorescence and absorption measurements on an X-ray Absorption Spectroscopy beamline[END_REF]. Successive isothermal plateaus corresponding to the green dots on Figure 6 will be performed under different atmospheres.

To fulfill by the safety rules, a 1 mg sample (1x1x0.1 mm 3 ) with a maximum activity of 14 Bq will be prepared and characterized through SEM before the experiment.

Additional thermal treatments

Additional thermal treatments will be performed in the DURANCE loop located at the UO2 laboratory, CEA Cadarache, France (green dots on Figure 6). During these tests, the sample is placed in a metallic Mo crucible and heated in an induction furnace where the atmosphere is controlled by a gas mixture (H2O/H2). The temperature is measured and controlled by a thermocouple placed below the crucible and the oxygen partial pressure is monitored thanks to a MicroPoas probe (GENAIR) [START_REF] Pontillon | Fission gas release under normal and off-normal conditions: New analytical device implemented at the CEA Cadarache[END_REF].

The samples will be then characterized through XAS, SEM-EDX and optical microscopy to complete the results obtained by in-situ XAS experiments.

Conclusion

An innovative method has been developed in order to study volatile FP speciation under SA conditions.

SPS syntheses

Simplified SIMFUEL samples containing Cs, Mo, Ba and Zr have been synthetized through SPS (1200°C, 5 min). Preliminary characterizations indicated that Cs is indeed present inside the samples after sintering. Further analyses, notably ICP-MS and SEM-EDX on polished samples, have to be performed in order to determine accurately the amount and initial chemical state of Cs.

The final density of the samples was roughly estimated to vary from 68 to 90 %dth depending on the mix. This is consistent with the observed microstructure confirming the sintering was uncomplete. Nevertheless, measurements of the density by the Archimedes method have to be performed to get more accurate data.

Thermodynamic calculations

Thermal treatments conditions, representative of a SA, were determined thanks to thermodynamic calculations using FACTSAGE coupled with the TAF-ID.

In-temperature XAS

Predominance diagrams indicating the FP chemical states in the different conditions will be probed by intemperature XAS analyses under controlled atmosphere performed on BM16 beamline at the ESRF. A specific sample holder and glovebox have been designed to handle and confine the sample during the experiment. This kind of analyses will enable the observation of potential intermediate compounds formed at high temperatures and improve the current knowledge concerning FP speciation.

Additional thermal treatments

Additional thermal treatments will be performed in DURANCE in order to complete the in-situ characterizations. The evolution of the samples' microstructure and the distribution of FP in the matrix will notably be studied by SEM-EDX.

Figure 1 :Figure 2 :

 12 Figure 1: SEM-SE image of the pre-reduced UO2 powder (mean ECD 30 µm)

Figure 3 :

 3 Figure 3: SEM-SE image of Cs nitrate and uranate after 24h at 700°C under air (mean agglomerate size 194.1 µm)

  5 kN (~88 MPa) and fill the chamber with Ar • Heat to 1200°C at 200°C/min • Hold at 1200°C for 5 min • Cool down at 200°C/min • Release the pressure The pellets were recovered and roughly polished to remove the graphite deposit on the surface. Their density was approximated by measuring their dimensions and compared to the theoretical density of the mix (equation (2)). The values varied between 68 and 90 %d 𝑡𝑡ℎ depending on the mix considered: samples of batch n°4 have a density of 86.1 ± 4.3 %d 𝑡𝑡ℎ . d 𝑡𝑡ℎ = ∑ 𝑑𝑑 𝑖𝑖 𝑛𝑛 𝑖𝑖 𝑤𝑤𝑤𝑤% 𝑖𝑖 is the theoretical density of the mix, d 𝑖𝑖 and wt% 𝑖𝑖 are the theoretical density and the weight percent of the chemical i in the mix, m, h and r are the mass, height and radius of the pellet resp.

  European Review Meeting on Severe Accident Research -ERMSAR-2017 Warsaw, Poland, 16-18 May 2017

Figure 5 :

 5 Figure 5: SEM-SE image of fractured surface from batch n°1. A melted area containing only Mo and is clearly visible.

  European Review Meeting on Severe Accident Research -ERMSAR-2017 Warsaw, Poland, 16-18 May 2017 Source Term issues and the Paper N° (order as in the announcement)

Figure 6 :

 6 Figure 6: Mo-Cs-U-O predominance diagram, FACTSAGE coupled with the TAF-ID

Figure 7 :

 7 Figure 7: In-temperature XAS apparatus with a heated wire as a sample holder (left) and newly designed sample holder (right): the sample is maintained on the wire by two crimped nuts. A thermocouple is placed as close as possible of the sample. A hole in the nuts enables the X-ray beam to interact with the sample.The experiment will be performed in High Energy Resolution Fluorescence Detection (HERFD) on FAME-UHD beamline (BM16) at the ESRF. XANES spectra will be recorded at the Cs LIII edge (5.012 keV). A Crystal

Table 1 : Composition of the different batches synthetized (concentrations in wt%)

 1

	Batch	UO2	Cs2MoO4 Cs2UO4	BaZrO3 BaMoO4	MoO3
	1	99.1	0.9
	2	98.4	0.9	0.7
	3	99.8	1.2
	4	98.1	1.2	0.7
	5	98.5	1.2		0.3

Acknowledgement

Authors are indebted to the GENTLE program for funding the syntheses of the samples at the JRC as well as to the JRC-Karlsruhe staff for their support and valuable contribution to this work. We greatly thank M. Ernstberger for the SEM-EDX results as well as P. Valenza for his help with the modifications of the intemperature XAS device. Finally, we acknowledge P. Colomp for his contribution to the safety analysis at the ESRF.