E A Hodille

Rémi Delaporte-Mathurin

J Denis

M Pecovnik

E Bernard

Y Ferro

R Sakamoto

Y Charles

J Mougenot

A De Backer

C S Becquart

S Markelj

C Grisolia

Modelling of hydrogen isotopes trapping, diffusion and permeation in divertor monoblocks under ITER-like conditions

Keywords: Hydrogen, materials, plasma material interactions, modelling Modelling of hydrogen isotopes trapping, diffusion and permeation in divertor monoblocks under ITER-like conditions2

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

Introduction

In a tokamak, the plasma facing components are exposed to intense fluxes of particles (≈ 10 24 m -2 s -1 ) and heat (≈ 10 MWm -2 ) [START_REF] Pitts | Physics basis for the first iter tungsten divertor[END_REF]. The divertor will experience the most intense plasma wall interaction. To sustain such harsh environment in steady state condition (minutes long plasma discharges), a tungsten (W) armor has been chosen which is actively cooled by a water navigating at 343 K in a cooper alloys (CuCrZr) tube.

In order to reduce mechanical incompatibilities between both materials [START_REF] Bui | The thorn singularity of a bimaterial solid in thermoelasto-plasticity[END_REF], a Cu inter layer is used. This is the design of the divertor plasma facing components of WEST [START_REF] Missirlian | The west project: Current status of the iter-like tungsten divertor[END_REF] and ITER.

This intense particle flux leads to the implantation and further trapping of hydrogen isotopes in the divertor monoblocks increasing the content of hydrogen in the wall elements of the tokamak. This causes two mains issues: i) A fusion plasma contains ideally an equal mixture of tritium and deuterium. The tritium being a radioactive element, the safety authority limits its overall inventory in the vacuum chamber of ITER to 700 g. ii) The fuel retained in the wall can act as an unpredictable source of plasma particles that can affect the core plasma: an uncontrolled outgassing from the wall element, especially during long pulses discharges (>minutes) can cause uncontrolled increase of the plasma density leading ultimately to its termination [START_REF] Grisolia | Plasma wall interaction during long pulse operation in tore supra[END_REF]. It is thus crucial to understand the mechanisms leading to fuel retention, release and permeation in plasma facing materials.

This work tackles this issue using macroscopic diffusion/reaction models based on the McNabb and Foster equations [START_REF] Mcnabb | A new analysis of the diffusion of hydrogen in rion and ferritic steels[END_REF]. In this particular study, the code MHIMS [START_REF] Hodille | Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials[END_REF][START_REF] Hodille | Simulations of atomic deuterium exposure in self-damaged tungsten[END_REF][START_REF] Hodille | Study of hydrogen isotopes behavior in tungsten by a multi trapping macroscopic rate equation model[END_REF][START_REF] Pečovnik | New rate equation model to describe the stabilization of displacement damage by hydrogen atoms during ion irradiation in tungsten[END_REF] (Migration of Hydrogen Isotopes in Materials) and the code FESTIM [START_REF] Delaporte-Mathurin | Finite element analysis of hydrogen retention in iter plasma facing components using festim[END_REF][START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF][START_REF] Delaporte-Mathurin | Influence of interface conditions on hydrogen transport studies[END_REF] (Finite Element Simulation of Tritium in Materials) are used. The MHIMS code is used to render the dynamic of hydrogen isotopes retention in one dimension (1D) from the subsurface up to the bulk during both pulsed and continuous plasma. The FESTIM code takes inputs from the MHIMS simulations and is used in a two dimension (2D) geometry to render the behaviour of hydrogen isotopes during continuous plasmas in a complex multi-material components.

Plasma facing material model: W divertor target

Brief model description

The model that describes HI migration in the plasma facing materials is based on the McNabb and Foster equations [START_REF] Mcnabb | A new analysis of the diffusion of hydrogen in rion and ferritic steels[END_REF]. In this model, the HI content is partitioned into two types of atoms: the mobile (interstitial) HI with the concentration c m (m -3 ) and the trapped HI with the concentration c t,i (m -3 ) for the trapping site of type i. In the classical McNabb and Foster model, the equations governing the time and spatial variation of concentration of both types of particles are recalled here:

∂c m ∂t = ∇ • D(T ) ∇c m + S - i ∂c t,i ∂t (1)
∂c t,i ∂t = ν t,i (T )c m (n i -c t,i ) -ν dt,i (T )c t,i (2)
where S (m -3 s -1 ) is a volumetric source of HI (from the implantation). The other parametres are defined further in this section. The local HI concentration (m -3 ) is calculated as c HI = c m + i c t,i . A Fickian diffusion is assumed for the mobile particles (first term on the right hand side of equation 1), which is characterised by the diffusion coefficient of HI in the material D(T ) (m 2 s -1 ). The ones used in this study are reported in [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF] for W, Cu and CuCrZr. Since the simulated plasma is composed of D atoms, the diffusion coefficient in each material is scaled by 1/ √ 2 [START_REF] Vineyard | Frequency factors and isotope effects in solid state rate processes[END_REF]. The mobile interstitial HI atoms can interact with different types of defects and become trapped HI atoms. Each defect kind i has a concentration n i (m -3 ). According to several atomic-scale simulations, defects such as monovacancies [START_REF] Fernandez | Hydrogen diffusion and vacancies formation in tungsten: Density functional theory calculations and statistical models[END_REF][START_REF] Heinola | Hydrogen interaction with point defects in tungsten[END_REF] or vacancy clusters [START_REF] Hou | Predictive model of hydrogen trapping and bubbling in nanovoids in bcc metals[END_REF] can retain several HIs with binding energies depending on the quantity of HI in the defects. These calculations led to a modification of the equation 2 of the classical model to take into account this multi-trapping aspect [START_REF] Schmid | Transport of hydrogen in metals with occupancy dependent trap energies[END_REF][START_REF] Hodille | Study of hydrogen isotopes behavior in tungsten by a multi trapping macroscopic rate equation model[END_REF].

The interactions between HI and trapping sites are described by two processes: detrapping and trapping. The detrapping rate constant ν dt,i (T ) (s -1 ) (second term of the right hand side of equation 2) is given by:

ν dt,i (T ) = ν 0,dt,i exp - E dt,i k B T (3)
with k B = 8.617×10 -5 eV/K, the Boltzmann constant, T (K) the temperature of the material, ν 0,dt,i = 10 13

s -1 the detrapping pre-exponential factor [START_REF] Fernandez | Hydrogen diffusion and vacancies formation in tungsten: Density functional theory calculations and statistical models[END_REF] and E dt the detrapping energy (eV) which depends on the amount of H in the defect (usually, it is calculated that the detrapping energy decreases with the number of trapped H). The trapping rate constant ν t (T ) (m 3 s -1 ) (first term of the right hand side of equation 2) is given by:

ν t,i (T ) = ν 0,t,i exp - E t k B T (4)
with ν 0,t,i (m 3 s -1 ) the trapping pre-exponential factor and E t the trapping energy. Assuming immobile defects and trapping energies equal to the diffusion energy, one can obtain [START_REF] Waite | Theoretical treatment of the kinetics of diffusion-limited reactions[END_REF]:

ν t (T ) = 4π(R t,i + λ)D(T ) (5)
with D(T ) (m 2 s -1 ) the diffusion coefficient of HI atoms in the lattice, λ (m) the capture radius of H in interstitial sites and R t,i the trapping radius of the trapping site i.

Model geometry

In WEST and ITER, the divertor targets are made of plasma facing units constituted of monoblocks composed of three materials:

(i) W as plasma facing material.

(ii) Cu as interlayer material.

(iii) CuCrZr alloy as material for the cooling pipes.

The geometry of such monoblock is shown in [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF] and in figure 1. Its thickness is e = 12 mm and its width is w = 28 mm. On the one hand, the ions from the plasma are implanted in the first nanometres in which specific trapping processes can occurs such as the formation of the super-saturated layer (SSL) [START_REF] Gao | Deuterium supersaturation in low-energy plasma-loaded tungsten surfaces[END_REF][START_REF] Gao | Hydrogen atom-ion synergy in surface lattice modification at sub-threshold energy[END_REF]. On the other hand, the retention in such a complex component requires a 2D geometry especially when dealing with high surface temperatures [START_REF] Delaporte-Mathurin | Finite element analysis of hydrogen retention in iter plasma facing components using festim[END_REF]. The MHIMS code which includes surface and subsurface models is well suited for simulations at the nanometremicrometre scale, while the FESTIM code is well suited for simlulations of transport in multidimensional and multi-material systems at the scale of the entire component. Thus, two sets of simulations will be performed: (i) 1D MHIMS simulations in the W part and (ii) 2D FESTIM simulations with the complete geometry of figure 1.

The MHIMS simulations are used to determine the amount of hydrogen isotopes (HI) retained in the supersaturated layers and to compare the HI retention during cycled and continuous plasma exposures. In these cases, only the W armour is simulated in 1D along the axis of symmetry of the monoblock: in that direction, the W width is 7.2 mm.

The FESTIM simulations are used to determine the total retention in the monoblocks, the permeation flux to the coolant and the desorption from the toroidal edges of the monoblocks. The 2D geometry of the monoblock presenting an axis of symmetry, only half of the monoblock is modelled. 1 and correspond to an attached plasma scenario according to SOLPS simulations [START_REF] Pitts | Physics basis for the first iter tungsten divertor[END_REF]. We consider a pure deuterium plasma so that the incident particles are deuterium atoms/ions. Table 1 gives the incident flux of particles but the source of implanted particles has to be corrected by the reflection coefficient, r ion . This reflection is estimated with the binary collision approximation (BCA) code SRIM (Stopping Range of Ions in Matter) [START_REF] Ziegler | The Stopping and Range of Ions in Matter[END_REF]. For deuterium ion impacting a W surface with 45 • with respect to the normal to the surface, the reflection coefficient is calculated to be ≈ 0.64.

The highest energy is 115 eV, which is below the threshold for physical sputtering of W [START_REF] Sugiyama | Sputtering of iron, chromium and tungsten by energetic deuterium ion bombardment[END_REF]. Thus, sputtering and redeposition of W surfaces atom are not accounted for in the current model. However, the kinetic energy distribution of the ions impacting the surfaces may not be mono-energetic and ions with kinetic energy above the sputtering threshold may interact with the surfaces of the plasma facing components, leading to sputtering and redeposition. The impact of sputtering and redeposition on the fuel retention estimation Table 1: Exposure conditions for the different scenarios considered during the steady-state plasma (SSP):

φ SSP heat (MWm -2 ),φ SSP ion (m -2 s -1
), E SSP ion (eV), the surface temperature T SSP surf (K) and the temperature at the rear of the W armour T SSP rear . The surface temperature T SSP surf resulting of such exposition is calculated using the formula in [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF].

scenario φ SSP heat φ SSP ion E SSP ions T SSP surf T SSP rear I 10 MWm -2 5.0×10 23 m -2 s -1 115 eV 1456 K 552 K II 5 MWm -2 3.5×10 23 m -2 s -1 75 eV 870 K 436 K III 1 MWm -2 1.0×10 23 m -2 s -1 50 eV 435 K 347 K
is discussed in more details in section 4.2. Also, it could be possible that the W surface inside the vessel is slightly covered with a native oxide layer, source of trapping sites. Such a oxide layer is not modelled here as one can consider it is sputtered by the D ions, even at 50 eV. The surface temperature is calculated using the equation ( 8) in [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF]: T surf = 1.1 × 10 -4 φ heat + 323 K with φ heat the heat flux (Wm -2 ). The temperature field is computed using the heat transport parameters that can also be found in [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF].

In the 1D-MHIMS simulations, two types of exposure are simulated:

-Cycled plasma: each cycle lasts 1,600 s with a 100 s plasma ramp-up, a 400 s steady-state plasma (SSP), a 100 s plasma ramp-down and a 1,000 s waiting phase before the next cycle. 25 cycles are simulated so that the duration of steadystate plasmas corresponds to 10,000 seconds.

During the ramp-up/down, the incident energy, implantation flux and temperatures (table 1) evolve linearly up to/from their steady-state values from/to its initial values (figure 2). -Continuous plasma: the exposure conditions in table 1 are kept constant for 10,000 s.

During the plasma exposure of the different scenarios, the temperatures at the rear of the 7.2-mm-thick W armour are taken from the FESTIM thermal module and are 552 K, 436 K and 347 K for the scenario I, II and III respectively. In the 2D-FESTIM simulations, only a 10,000 s continuous plasma is simulated. The implanted particles coming from the plasma act as a volumetric source of mobile particles deposited in the implantation zone. It can be simulated by S in the right hand side of equation 1 with S = f (x)(1r ion )φ ion . The implantation profile f (x) (in m -1 ) gives the distribution of ions/atoms once they have lost all their kinetic energy to the target through collisions to the atoms (nuclear stopping power) or friction with the electronic cloud (electronic stopping power). At that point, the particles start diffusing with a Fickian diffusion, which is handled by the current model. Thus, f (x) is calculated using BCA methods as for the reflection coefficient, which can take into account the various energy dissipation processes (nuclear and electronic stopping power) that leads to stopping of ions in the metal. For the incident energies used here (table 1), the implantation zones f (x) do not extent further than 10 nm below the surface. In 1D MHIMS simulations, it is possible to handle such volumetric source and the HI transport in the millimetre scale using a 500-cells mesh refined at the plasma exposed side. However, such a refinement in a 2D geometry is more problematic if one wants to stay in reasonable computation time. Thus, we decided to first simulate the exposure with the 1D MHIMS simulations and a volumetric source. Then, the results (see section 2.5) of these simulations are provided as input of 2D FES-TIM simulations to determine the full D concentration field in the 2D. For scenarios I and II for which the temperatures are high, the deep migration at the millimetre-scale requires the use of the full geometry of the monoblock (half considering the symmetry) with a mesh of ≈ 60,000-elements. For scenario III, the HI migration does not extend deep in the bulk and only the top 0.5 millimetre are simulated as the simulation of the rest of the component is not necessary (see section 3).

Trapping parameters

2.4.1. in the subsurface: supersaturated layer It has been shown experimentally that when HIs are implanted in W with high fluxes at low temperature, even below the displacement threshold energy, the formation of a so-called super-saturated layer (SSL) can be observed. In the SSL, the concentration of H atoms can reach about 0.01 -0.1 atomic fraction (at.fr.) in the implantation zone (up to 10 nm) [START_REF] Kh | Depth distribution of deuterium in single-and polycrystalline tungsten up to depths of several micrometers[END_REF][START_REF] Gao | Deuterium supersaturation in low-energy plasma-loaded tungsten surfaces[END_REF][START_REF] Gao | Hydrogen atom-ion synergy in surface lattice modification at sub-threshold energy[END_REF]. The formation of the SSL is theoretically understood as the formation of mono-vacancies by the introduction of a large quantity of hydrogen in interstitial position, while it is energetically more favourable for the implanted H to sit in vacancies [START_REF] Hodille | Hydrogen super-saturated layers in h/d plasma loaded tungsten: A global model based on thermodynamics, kinetics and density functional theory data[END_REF]. This SSL is only formed in the implantation zone of the ions where the ions transfer their energy to the lattice. The explanation proposed in [START_REF] Hodille | Hydrogen super-saturated layers in h/d plasma loaded tungsten: A global model based on thermodynamics, kinetics and density functional theory data[END_REF] and [START_REF] Gao | Hydrogen atom-ion synergy in surface lattice modification at sub-threshold energy[END_REF] is that the implantation zone is the only place where the kinetic barriers can be overcome to reach the thermodynamic equilibrium. According to the model described in [START_REF] Hodille | Hydrogen super-saturated layers in h/d plasma loaded tungsten: A global model based on thermodynamics, kinetics and density functional theory data[END_REF], the SSL formation is conditioned upon the exposure conditions (temperature and incident flux of hydrogen). According to this model, only scenario III would lead to the growth and stabilisation of the SSL (low temperature and high flux). For scenarios I and II, the temperature is too high. However, one can consider that the SSL can still be present either thanks to the previous plasma pulses or thanks to the ramp-up/down phase during which the conditions to form the SSL are met.

The trapping parameters, i.e. the detrapping energies, E dt of HI from the mono-vacancies and the trapping radius R t are given in table 2. Several DFT calculations are available for the binding of HI atoms with a mono-vacancy in W [START_REF] Fernandez | Hydrogen diffusion and vacancies formation in tungsten: Density functional theory calculations and statistical models[END_REF][START_REF] Heinola | Hydrogen interaction with point defects in tungsten[END_REF][START_REF] Hou | Predictive model of hydrogen trapping and bubbling in nanovoids in bcc metals[END_REF][START_REF] You | Dissolving, trapping and detrapping mechanisms of hydrogen in bcc and fcc transition metals[END_REF]. All these studies agree on the fact that the detrapping of HI from a mono-vacancy ranges from about 0.8 eV for six trapped HI atoms to 1.5 eV for one trapped HI atom. In this work, the detrapping energies calculated by Fernandez et al. [START_REF] Fernandez | Hydrogen diffusion and vacancies formation in tungsten: Density functional theory calculations and statistical models[END_REF] are chosen to described the energetic of HI trapping in mono-vacancies.

In the current model, mono-vacancies are assumed to have constant density with time, restricted to the implantation zone (i.e, up to 10 nm). Their concentration in the SSL is set to 5×10 -3 at.fr. (3.2×10 26 m -3 ). Thus, once filled with 6 HI, the HI concentration in the SSL reaches ≈ 3×10 -2 at.fr. (1.9×10 27 m -3 ) as observed experimentally.

Table 2: Trapping parameters for the SSL: detrapping energies E dt (eV) and capture radius R t (nm). The detrapping energies are of HI from mono-vacancies are directly taken from [START_REF] Fernandez | Hydrogen diffusion and vacancies formation in tungsten: Density functional theory calculations and statistical models[END_REF]. The trapping radius is calculated as:

R t (V m ) = a √ 3 4 + 3m 8π 1 3 -3 8π 1 3 [26]. with a = 316 pm, the W lattice constant. V 1 D n E dt R t V 1 D 1 1.43 eV V 1 D 2 1.42 eV V 1 D 3 1.25 eV 0.137 nm V 1 D 4 1.17 eV V 1 D 5 1.11 eV V 1 D 6 0.86 eV 2.4.2.
In the bulk In the W armour, in addition to the SSL traps, we add two intrinsic traps with detrapping energies of 0.85 eV and 1.00 eV which have proven to be good trapping parameters to describe the bulk defects of polycrystalline W [START_REF] Hodille | Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials[END_REF]. A detrapping energy of 0.85 eV can be attributed to dislocation lines [START_REF] De Backer | Hydrogen accumulation around dislocation loops and edge dislocations: from atomistic to mesoscopic scales in BCC tungsten[END_REF] or metallic impurities like Fe, Ni or Cu [START_REF] Lu | A review of modelling and simulation of hydrogen behaviour in tungsten at different scales[END_REF]. A detrapping energy of 1.00 eV can be attributed to grain boundaries (GB) [START_REF] Zhou | Investigating behaviours of hydrogen in a tungsten grain boundary by first principles: from dissolution and diffusion to a trapping mechanism[END_REF][START_REF] Bakaev | Ab initio study of the stability of h-he clusters at lattice defects in tungsten[END_REF]. The GB used in the calculations of these references are symmetrical because they were obtained with periodic DFT codes. A more realistic energetic representation of GBs would be a continuous distribution function of activation energies as mentioned in [START_REF] Piaggi | Hydrogen diffusion and trapping in nanocrystalline tungsten[END_REF], where H trapping at random GB in nanocrystalline W is studied with molecular dynamics (with a bond order potential [START_REF] Li | Modified analytical interatomic potential for a w-h system with defects[END_REF]). Such type of continuous distribution of detrapping energies is not yet implemented in the model and we used instead the detrapping energy that allows reproducing experimental TDS measurements. The trapping parameters for these intrinsic traps are reported in 3. Their capture radii are set to be λ and there concentration to 10 -4 at.fr.. For the 2D FESTIM simulations, HI may reach the cooling system and permeate through W, Cu and CuCrZr. The transport and trapping parameters for Cu and CuCrZr used for these simulations are the same as the one given in [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF]. At material interfaces, the conservation of the chemical potential is guaranteed with the model described in [START_REF] Delaporte-Mathurin | Influence of interface conditions on hydrogen transport studies[END_REF]. This interface model required the solubility of HIs in the different materials taken from this same reference.

Boundary conditions

For the 1D MHIMS simulations, the boundary conditions (BC) for the mobile particle concentration Table 3: Trapping parameters for the intrinsic traps: detrapping energies E dt (eV), trap concentration (n i ) and capture radius R t (nm). The trap concentrations are expressed given in m -3 and also in atomic fraction, i.e. normalised to the tungsten concentration: ρ W = 6.338 × 10 28 m -3 .

E dt

n i R t defect type 0.85 eV 6.338 × 10 24 m -3 (10 -4 at.fr.) λ Dislocations, Metallic impurities [27, 28] 1.00 eV 6.338 × 10 24 m -3 (10 -4 at.fr.) λ Grain boundaries (GB) [START_REF] Zhou | Investigating behaviours of hydrogen in a tungsten grain boundary by first principles: from dissolution and diffusion to a trapping mechanism[END_REF] are: (i) On the plasma facing surface, a Dirichlet condition: c m (0, t) = 0.

It represents an instantaneous desorption of all particles arriving at the surface. (ii) On the rear side of the W armour, a Dirichlet condition: c m (7.2 mm, t) = 0. The BC for the rear side of the W armour is the closest conditions one can use with MHIMS to simulate a realistic permeation behaviour of HIs in the W armour.

For the 2D FESTIM simulations in the semimonoblock, the BC for the mobile particle concentration are:

(i) On the plasma facing surface (Γ top , a Dirichlet condition:

c m = c MHIMS m (X imp ).
With X imp being the implantation depth of the HI atoms. (ii) On the surface of the cooling pipe (Γ coolant , a condition based on a recombination coefficient:

-D CuCrZr ∇c m • n = K rec c 2 m with D CuCrZr (m 2 s -1
) the diffusion coefficient of HI in CuCrZr, K rec (m 4 s -1 ) the recombination coefficient taken from [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF] and n the normal to the surface. (iii) On the right side of the semi-monoblock, a homogeneous Neumann condition: ∇c m • n = 0 to account for symmetry conditions. (iv) On the left side of the semi-monoblock (Γ lateral ), a Dirichlet condition: c m = 0 that simulates the desorption of HI from the toroidal edges of the monoblock, in the toroidal gap. For a given scenario, first, MHIMS calculates the concentration of mobile HI in the implantation zone, c MHIMS m (X imp ) with the ions coming from the plasma being treated as a volumetric source for the different scenario.

Then, FESTIM uses this calculated concentration as an input to calculate the D concentration field in the 2D geometry of the monoblock.

Simulations Results

1D simulations

First, we present the results of the MHIMS simulations especially the impact of simulating cycled (with the evolution given in figure 2) or continuous plasma. There is a big difference between the depth profiles before the ramp down (black dashed lines) and at the end of the cycles for scenarios I and II. This is explained by the retention during the 100 seconds of plasma ramp-down during which the temperature decreases while the incident flux is not zero: trapping takes over detrapping increasing the filling ratio of the traps. Such ramp-down retention has already been briefly described in [START_REF] Delaporte-Mathurin | Finite element analysis of hydrogen retention in iter plasma facing components using festim[END_REF] and it will be investigated further in the following. It is also important to note that the quantity of D retained during this part of the plasma is independent of the cycles as one cannot distinguish the D profiles in this first mm.

The difference between the cycles is observed in the bulk part of the depth profile which shows the migration further and further in the bulk as the number of cycles increases. In addition, for scenarios I and II, the D concentration tends to increase close to the rear part of the W where the temperature is cooler than at the front exposed to the plasma.

Even though the retention in the SSL after the ramp down is similar for the three scenarios, it is different during the steady-state plasma (black dashed lines in figure 3). In Scenario III, the vacancies in the SSL are almost exclusively (99%) filled with 6 D during the steady state plasma. In scenario II, there are mostly VD 2 (42%), VD 3 (32%) and VD 4 (15%) with a few percents of empty vacancies. In scenario I, there are almost only empty vacancies (99%). Thus, during scenarios II and III, since most of the vacancies are filled with at least 1 D, the structure of the defects in this SSL should not change that much as the presence of D inhibits the migration [START_REF] Fernandez | Hydrogen diffusion and vacancies formation in tungsten: Density functional theory calculations and statistical models[END_REF] and thus coalescence of vacancies. For scenario II, the few percents of empty mono-vacancies could possibly move and aggregate with filled vacancies as shown in [START_REF] Heinola | Modelling of the effect of elms on fuel retention at the bulk w divertor of jet[END_REF]. For scenario I, almost all vacancies are empty and one could expect them to eventually be annealed following the various annealing steps reported experimentally in [START_REF] Debelle | First temperature stage evolution of irradiation-induced defects in tungsten studied by positron annihilation spectroscopy[END_REF][START_REF] Zibrov | Vacancy cluster growth and thermal recovery in hydrogen-irradiated tungsten[END_REF]. According to Zibrov et al. [START_REF] Zibrov | Vacancy cluster growth and thermal recovery in hydrogen-irradiated tungsten[END_REF], the SSL would actually disappear at 1456 K during the steady-state plasma but could be reconstructed during the ramp-down as the temperature decreases. 4 shows the evolution of the D retention for the various scenarios and for cycled (solid line) and continuous (dashed line) exposures as a function of the implanted fluence. The spikes appearing periodically on the retention for scenarios I and II are the retention occurring during the ramp down phase which is removed during the next ramp-up phase. When focusing only on the evolution of the retention between these spikes, i.e. during the steady-state plasma, no obvious differences can be observed between cycled and continuous plasma, at least in the time dynamics: the retention due to the migration of D deeper and deeper in the bulk is the same regardless if the plasmas are cycled or not. Thus, the study of D retention and permeation for long exposure times can be studied with either types of plasma exposures.

cycle vs continuous plasma: D retention Figure

The evolution of the retention has been interpolated using the curve fit function of the python package scipy [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF][START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF] with a power law. For scenario III, the exponent of the power law is 0.5 showing a retention limited by the migration in the bulk. And as the surface temperature increases, the exponent deviates from 0.5 to reach 0.66. This effect is explained by the lower temperature at the rear side of the W armour which tends to increase the concentration in this area. In [START_REF] Hodille | Estimation of the tritium retention in ITER tungsten divertor target using macroscopic rate equations simulations[END_REF], we already showed an evolution of the exponent but it was going below 0.5. We were looking at the evolution of the retention at the end of the cycle which accounts for the retention during the ramp-down which is not taken into account in figure 4. This is because we were not taking the evolution of the retention only during the steady-state plasma but we were accounting for the retention during the ramp-down. Doing so with the present data would also make the retention evolve with a power law below 0.5. However, if one wants the exponent to make sense with the bulk migration dynamic, one should focus only on the retention during the steady-state plasma.

3.1.3. Ramp-down retention/Ramp-up release Figure 5 shows the evolution of the D retention during the 10 th cycle of scenario II. During the ramp down phase (first 100 seconds), the retention decreases drastically due to the temperature increase. This leads to a release during the ramp down highlighted by the arrow RU rel in figure 5. During the steady-state plasma, the retention keeps decreasing in the next 100 seconds before increasing again. When the plasma ramp down starts, the retention goes up during the ramp down which is highlighted by the arrow RD ret in figure 5. This increase of the D retention appears because the PFC is actively cooled: the particle flux and the temperature decreases with the same behaviour (figure 2), and eventually the temperature is low enough and the flux of particle is still high enough to favour trapping over detrapping. If the fluxes were immediately zero after the steady state plasma or if the temperature decreased much slower than the flux (inertial PFC for instance), no ramp down retention would be observed. Figure 6 reports the quantity of D released/retained during the ramp-up/down in the three scenarios. These quantities are computed as shown on figure 5 (see the arrows). No release during the ramp-up is reported for the first cycle as the material is empty at the beginning of the simulation: for any scenario, it would be negative (uptake of D) and would not add anything significant to the interpretation of the dynamic of the D retention.

For scenario III, the retention during the ramp down tends to decrease and the release during the ramp-up is negative meaning there is an uptake of D during this phase. This behaviour is expected for a cold scenario with small variations of temperature.

For scenarios I and II, the retention during the ramp down is positive as well as the release during the ramp-up. In addition, both are very close, the difference being that during the first 100 seconds of the steady-state plasma, the retention keeps decreasing (see figure 5). Thus, the retention during the ramp down phase, localised close to the exposed surface of the PFC, is removed during the ramp-up phase, as explained in the previous section. This dynamic retention stays constant of the order of 10 20 Dm -2 during all cycles and may act as a source of fuel for the plasma during the ramp down.

The WEST divertor is composed of 480 PFUs and assuming there are two monoblocks per PFU exposed to scenarios I and II (outer and inner), one can estimate that the dynamic retention during the rampup/down represents about 2 × 10 20 D atoms. Taking a WEST plasma of 15 m 3 and a plasma density of 9 × 10 19 m -3 (high power scenario in [START_REF] Bucalossi | The west project: Testing iter divertor high heat flux component technology in a steady state tokamak environment[END_REF]), the retention/release during the ramp down/up represents about 15% of the total amount of D in the plasma. The ITER divertor is composed of 54 cassettes with 16 PFUs on the inner side and 22 PFUs on the outer sides. Assuming the is one monoblock per PFU exposed to scenarios I and II, the total quantity of D retained/released during the ramp down/up is of the order of 4×10 20 D atoms. Taking an ITER plasma of 830 m 3 with a density of 10 20 m -3 , the retention/release during the ramp down/up represents 0.5% of the total amount of D in the plasma. Thus, this retention and release could be a possible issue for the control of the source of fuel during the plasma ramp up especially in smaller tokamaks. Indeed, currently, there is no direct coupling between plasma transport codes such as SOLPS [START_REF] Wiesen | The new solps-iter code package[END_REF][START_REF] Kaveeva | SOLPS-ITER modelling of ITER edge plasma with drifts and currents[END_REF] or SOLEDGE3X-EIRENE [START_REF] Bufferand | Numerical modelling for divertor design of the west device with a focus on plasma-wall interactions[END_REF][START_REF] Bufferand | Threedimensional modelling of edge multi-component plasma taking into account realistic wall geometry[END_REF] and retention codes like the ones used here. Thus, no ionisation of the emitted particles is taken into account in the retention code, that would eventually increase the local plasma density and possibly affect the exposure conditions. However, First attempts to calculate recycling coefficient with retention codes have emerged [START_REF] Denis | Dynamic modelling of local fuel inventory and desorption in the whole tokamak vacuum vessel for autoconsistent plasma-wall interaction simulations[END_REF] with the goal of coupling SOLEDGE3X-EIRENE and the Desorption from Wall ElemEnts (D-WEE) module, which contains MHIMS to model the HI/material interactions.

It is important to note that the quantity given here are indicative and would likely change with the width of the heat/particle deposition on the divertor or with the duration of the ramp-up/down. For instance, a shorter ramp down may reduce the retention during the ramp down. On the contrary, a wider deposition of the fluxes on the divertor would increase the overall quantity of D retained during the ramp down.

2D simulations

The concentration of mobile particles, c MHIMS m (X imp ), calculated in the 1D MHIMS and given as input to FESTIM are reported in table 4. In the MHIMS simulations, the concentration of mobile particles reached 95% of the value reported in table 4 in less than a second. Thus, a constant concentration is chosen as input of the FESTIM simulations. One can see that the concentration of mobile D is higher for the coldest scenario. Indeed, the mobility of D is related to (i) the gradient of the concentration of mobile particles (i.e. to φ inc , and E inc ) and (ii) the diffusion coefficient (i.e. the temperature T surf ). This is calculated by the model with a volumetric implantation source (with MHIMS) and given as an input to FESTIM. However, one can estimate the concentration of mobile particles for a given set of exposure conditions assuming no (or negligible) gradient to the bulk and a full recycling of implanted D [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF][START_REF] Gao | Hydrogen atom-ion synergy in surface lattice modification at sub-threshold energy[END_REF][START_REF] Hodille | Hydrogen super-saturated layers in h/d plasma loaded tungsten: A global model based on thermodynamics, kinetics and density functional theory data[END_REF]:

c m = (1-rion)φincRimp D(T )
with R imp (m) the implantation depth of the ions. Since (i) from scenario I to III the diffusion coefficient D(T ) drops by 42, and (ii) the incident flux φ inc and the implantation depth R imp only decrease by 5 and 1.6 respectively, the above assumption on c m leads to a 5 time higher concentration of mobile particles for scenario III.

The conclusion of the previous sections concerning the retention, migration and permeation in the far bulk is that one can simulate either continuous or cycled plasma exposures without changing the results. The continuous plasma are much faster to run than cycled plasma with FESTIM so we decided to run 10,000 seconds continuous plasma exposures. First, the 2D D retention fields after 10,000 seconds of exposure are presented and then macroscopic data are extracted from the simulations.

3.2.1.

2D retention fields Figure 7 shows the deuterium retention in the half-monoblock after 10,000 seconds of plasma exposure for scenario I (a) and II (b). As already mentioned in [START_REF] Delaporte-Mathurin | Finite element analysis of hydrogen retention in iter plasma facing components using festim[END_REF][START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF][START_REF] Delaporte-Mathurin | Influence of interface conditions on hydrogen transport studies[END_REF], and as observed in the 1D depth profiles (figure 3), the presence of the thermal gradient implies that the maximum of D concentration is not at the exposed surface but near the cooling pipe. This is especially visible for scenario I (hottest scenario) for which the maximum of retention tends to follow the shape of the cooling pipes. The impact of the interface model is not seen on the D concentration of scenario II but it is seen for scenario I. The continuity of the chemical potential chosen here leads to a discontinuity of the concentration of mobile particles at the W/Cu and Cu/CuCrZr interfaces. Especially, as shown in [START_REF] Delaporte-Mathurin | Influence of interface conditions on hydrogen transport studies[END_REF], this tends to enhance the retention in the CuCrZr part.

According to the 1D simulations, the D migration depth does not exceed 500 µm for scenario III (figure 3(c)). Thus, for this scenario, simulating the D transport in the same geometry as for scenarios I and II is not necessary. Instead, we simulate a 500 µm×14 mm layer of W. One of the 14 mm boundary (Γ top ) is exposed to the plasma (Dirichlet BC see section 2.5) and D is allow to desorb instantaneously from one of the 500 µm boundary (Γ lateral ). Both are Dirichlet BC as explained in see section 2.5. The temperature in this 500 µm×14 mm rectangle is set to be the temperature of scenario III. The resulting D retention field after a 10,000 seconds exposure is shown in figure 8 (it does not show the full 14 mm width but it is zoomed on the corner where D is allowed to desorb). As in the 1D simulations, far from the desorbing edge, the D concentration is constant up to the migration depth which is about 300 µm. However, close to the toroidal edge, the migration depth decrease, which is also seen for scenario II (figure 7(b)). This is because D is allow to desorb instantaneously from the toroidal edge to the gap between PFUs forcing the D concentration at this edge to 0. These edges effects stay marginal if one only focuses on the D retention. However, the desorption from the toroidal edges is a source of fuel recycling from the wall to the plasma which is estimated in the next section. . The SSL is very close to the surface (10 nm) and the retention in it is not influenced (or only marginally) by the 2D shape of the monoblock. Thus, we used the linear combination of equation 6 to calculate the D retention.

This quantity is shown in figure 9 for the three scenarios during a continuous 10,000 seconds plasma The thickness of the simulated material is 500 µm. The width is the same as for the other scenarios (14 mm) but for sake of clarity, the figure is zoomed in the corner from which the D is allow to desorb (Γ lateral ). The temperature is 435 K. exposure in addition to the SSL contribution to the SSL. The SSL contribution does not evolve during the continuous exposure: the SSL is located in the first 10 nm below the surface and is almost instantaneously saturated. For scenario II and III, this contribution is > 10 15 D per monoblock (≈ 10 18-19 Dm -2 ) and represent few % of the total inventory at the beginning of the simulations but becomes rapidly negligible on the total inventory. It is however important to assess the content of D in the SSL as it is a close source of fuel and transient like ELMs could modify that content [START_REF] Heinola | Modelling of the effect of elms on fuel retention at the bulk w divertor of jet[END_REF] and perturb the recycling from the wall to the plasma [START_REF] Schmid | Diffusion-trapping modelling of hydrogen recycling in tungsten under ELM-like heat loads[END_REF]. For scenario I, the content of D in the SSL is completely negligible compared to the retention in the bulk which is expected considering the high temperature in this scenario. Despite the much lower migration depth, the inventory is maximum for the lowest temperature. This is consistent with previous simulation work in 1D [START_REF] Hodille | Estimation of the tritium retention in ITER tungsten divertor target using macroscopic rate equations simulations[END_REF][START_REF] Sang | Modeling of fuel retention in the upper tungsten divertor of east from attached to detached divertor plasma[END_REF] and 2D [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF] geometry. The evolution of the retention has also been fit with power laws. For scenario III, the retention evolves as a square root of time which agrees well with the steady-state model presented in [START_REF] Hodille | Simulations of atomic deuterium exposure in self-damaged tungsten[END_REF] for migration limited bulk retention. For scenario I and II, the power law slightly deviates from the square root of time as shown on figure 4. However, the exponent obtained from the 1D simulations is further away from the square root of time than the one obtained from the 2D simulations which is explained by the 2D complex migration of D around the cooling pipe, especially in scenario I (figure 7(a)).

3.3.2. Permeation flux to the coolant Figure 7(a) shows that D reached the CuCrZr pipes in 10,000 seconds. It leads to permeation of D from the plasma facing side to the cooling system creating a leak of D which is shown in figure 10(a). This permeation is observed only for the simulation of scenario I with the hottest temperature. After 10,000 seconds of continuous plasma, the quantity that is released to the water cooling system is calculated to be 1.35×10 14 D per monoblock with a flux ≈ 0.6 × 10 11 D/s per monoblock. This flux is still growing as the exposure time showing the same transient form as in [START_REF] Delaporte-Mathurin | Influence of interface conditions on hydrogen transport studies[END_REF] in which the permeation flux eventually stabilizes (at least for the 1D simulations). Thus, one could expect that the permeation flux would eventually reach a steady-state, though the retention is still growing after 10 7 seconds of plasma exposure in [START_REF] Delaporte-Mathurin | Influence of interface conditions on hydrogen transport studies[END_REF] making the time to reach the steady-state far beyond our simulation time scale.

In section 3.1.3 we assumed that one monoblock per PFU would be exposed to scenario I at the outer and inner targets: the total quantity released to the cooling system after 10,000 seconds is 2.8×10 17 D atoms and a permeation flux 1.2 × 10 14 D/s. If one extend these results to the ITER nuclear phase during which the plasma is a half deuterium/half tritium mixture, it means the permeation flux to the cooling system after 10,000 seconds of continuous exposure is estimated to 0.1 MBq/s (50 Bq/s per monoblock) for a total of 252 MBq after 10,000 seconds of continuous plasma. This would require tritium removal from the cooling system, if the monoblocks are exposed to such a high temperature for more than 10,000 seconds in total. It is important to note that the interface model assume perfection of the interfaces between both material: there is no traps which slows down the permeation processes. The impact of trap at the interface may be significant if the trap is energetically deep as shown in the appendix of [START_REF] Delaporte-Mathurin | Influence of interface conditions on hydrogen transport studies[END_REF] but it requires a parametrisation which is not yet available. The values given here is thus an upper limit for the permeation flux.

3.3.3.

Outgassing at the toroidal edge Figure 7 and figure 8 show that the D retention field is affected by the desorption from the toroidal edge of the monoblock Γ lateral . The desorption flux from these edges to the gap between monoblocks is calculated as 2e Γ lateral (-D(T )∇c m • n) dl and shown in figure 10(b) for the three scenarios. The growth seems to have two time constants: a fast one which contributes to the most of the desorption flux, and a slower one with a much smaller contribution. It means that most of the desorbed D is done at the upper corner of the toroidal edge (fast time constants) and the slow time constant is due to the slow migration in the bulk. The former process has a smaller contribution as D migrates toward colder area where D(T ) and thus the desorption flux is smaller (and the gradient of mobile particle is lower as well).

The formula for the desorption flux implies that the desorption flux is basically driven by the value of the gradient of c m , hence by the value of c MHIMS m (X imp ) (table 4) and by the value of the diffusion coefficient, hence the temperature. A first crude assumption would then means the desorbing flux is proportional to ∝ D(T surf )c MHIMS m (X imp ) which is validated by the ratio between the desorbing flux of the three scenario: the desorbing flux of scenario I is twice the desorbing flux of scenario II and 10 times the scenario of scenario III.

Considering that one monoblock per inner/outer PFU target is exposed to each of these scenarios, the total quantity of D desorbing from the toroidal edges is 1.55×10 17 D/s (5.8×10 -4 Pa.m 3 .s -1 ) for the WEST divertor and 3.32×10 17 D/s (1.2×10 -3 Pa.m 3 .s -1 ) for ITER. This amount is very low (2 order of magnitude below the implanted flux on the top surface of the monoblock) thus, the desorption from the toroidal edges is negligible with respect to the desorption from the front surface. We simulate a 2D geometry: no poloidal edges are simulated and it should be another source of desorption from the monoblocks, but still very moderate.

Discussion

Permeation short cut through the poloidal gaps

Zhou et al [START_REF] Zhou | A tritium permeation 'short cut' for plasma-facing components of fusion reactors[END_REF] showed experimentally that HI atoms can penetrate in the gaps between adjacent monoblock in the poloidal direction. In these gaps, HI are directly in contact with the CuCrZr cooling pipe, thus offering a short cut for HI to reach the cooling system. An estimate of tens of seconds was given in [START_REF] Zhou | A tritium permeation 'short cut' for plasma-facing components of fusion reactors[END_REF], which is much faster than the thousands of seconds required for HI to diffuse across the full monoblock as shown in figure 10 (a). However, a proper estimate of the quantity of HI reaching the coolant via that channel would require the knowledge of the gas pressure and/or the plasma density of HI in these gaps. Such estimate is not easy to do since HIs come from two different sources: (i) from the plasma penetrating the gaps, and (ii) from the HIs desorbing from the poloidal edges of the monoblocks. In principle, the former (i) can be determine by plasma transport and gas kinetic models, while the latter (ii) required 3D FESTIM simulations of H transport and desorption in the monoblocks. In addition, 3D simulations of the monoblock plus cooling pipe are also required to calculate the permeation flux to the cooling pipe coming from both the gap short cut and the monoblock. Due to the complexity of such simulations, and given the uncertainty of the HIs flux via the short cut, we here focus on the best case scenario regarding the HI permeation time scale to the coolant and we leave the above mentioned simulations for another complete work.

Impact of the sputtering/redeposition on the HI retention

The highest energy considered is 115 eV which is below the threshold for physical sputtering of W according to [START_REF] Sugiyama | Sputtering of iron, chromium and tungsten by energetic deuterium ion bombardment[END_REF]. However, in an edge plasma, the particles are not mono-energetic: there is a distribution given by the ionic temperature. Thus, the incident energies of the particles hitting the wall are given by this distribution affected by the presence of the magnetic sheath (that tends to accelerate the ions). One can expect to have some fast ions that may sputter the W armor, though at a very low rate (3×10 -4 atoms per incoming ions according to [START_REF] Sugiyama | Sputtering of iron, chromium and tungsten by energetic deuterium ion bombardment[END_REF]). To take into account this effect, one could use a similar workflow as described in [START_REF] Lasa | Multi-physics modeling of the long-term evolution of helium plasma exposed surfaces[END_REF] that gives a ionic temperature to a plasma sheath code to obtain the distribution of the ion energy on the wall. Then, a binary collision approximation (BCA) model and an impurity transport model are used to estimate the erosion and redeposition rates. The results presented in [START_REF] Lasa | Multi-physics modeling of the long-term evolution of helium plasma exposed surfaces[END_REF] suggest the presence of an erosion zone around the strike point (due to the He ions produced by the fusion reactions) and some of these sputtered W atoms are redeposited toward the private flux region. This erosion/redeposition may affect the retention calculated in this study, especially in the sub-surface region and the SSL.

Scenario I and II would most likely be representative of erosion zones (being at or close to the strike point). The exposure conditions of these scenarios do not allow the formation of the SSL [START_REF] Hodille | Hydrogen super-saturated layers in h/d plasma loaded tungsten: A global model based on thermodynamics, kinetics and density functional theory data[END_REF] (too high temperature) but we still assumed such SSL can be present due to past exposure (ramp up/down phases or colder irradiation conditions). Thus, if these scenarios are in an erosion zone, the SSL would be progressively removed without any possibility to be rebuilt except during ramp up/down phases. Consequently, the contribution of the SSL to the total inventory (dashed line of figure 9) would disappear, knowing it is already low compared to bulk HI inventory. For the plasma background calculated in [START_REF] Lasa | Multi-physics modeling of the long-term evolution of helium plasma exposed surfaces[END_REF], the net erosion flux is estimated to be of the order of 10 18 m -2 s -1 (0.1 monolayer per second). Taking this value as a first estimate and the thickness of 1 monolayer to be 1 Å, a plasma of 10,000 s would remove about 100 nm of material, removing completely the SSL.

Scenario III would most likely be representative of a deposition zone. For this scenario, the formation of the SSL is possible according to [START_REF] Hodille | Hydrogen super-saturated layers in h/d plasma loaded tungsten: A global model based on thermodynamics, kinetics and density functional theory data[END_REF] (high flux and low temperature). The SSL is formed in the implantation zone where the energy is deposited: the SSL will only grow in the first nanometers and the SSL will build up as the material is deposited. The final quantity of HI in the SSL will be linked to a tradeoff between the SSL traps (mono-vacancies) creation rate and the material deposition rate. If the material is very rapidely deposited compared to the creation of the SSL traps, the SSL will be thick but with a low amount of HI while if the creation of the SSL traps is fast compared to the material deposition, the SSL will be thin but with a high concentration of HI.

Defects formation and evolution

4.3.1. Super-saturated layer We modelled a SSL in the sub-surface by a fixed concentration of monovacancies up to 10 nm. However, such a SSL is induced by D ion irradiations [START_REF] Gao | Deuterium supersaturation in low-energy plasma-loaded tungsten surfaces[END_REF][START_REF] Gao | Hydrogen atom-ion synergy in surface lattice modification at sub-threshold energy[END_REF] and one should include the dynamics of defect formation in the model using trap creation equation as proposed in [START_REF] Ogorodnikova | Deuterium retention in tungsten in dependence of the surface conditions[END_REF][START_REF] Ogorodnikova | Ion-driven deuterium retention in tungsten[END_REF]. Such model exists in MHIMS [START_REF] Hodille | Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials[END_REF][START_REF] Hodille | Stabilization of defects by the presence of hydrogen in tungsten: simultaneous w-ion damaging and d-atom exposure[END_REF] and have even been upgraded to account for the stabilisation of traps by the presence of trapped H [START_REF] Pečovnik | New rate equation model to describe the stabilization of displacement damage by hydrogen atoms during ion irradiation in tungsten[END_REF][START_REF] Pečovnik | Experiments and modelling of multiple sequential MeV ion irradiations and deuterium exposures in tungsten[END_REF]. So far, the domain in which the SSL could be formed has been studied with thermodynamic models [START_REF] Fernandez | Hydrogen diffusion and vacancies formation in tungsten: Density functional theory calculations and statistical models[END_REF][START_REF] Hodille | Hydrogen super-saturated layers in h/d plasma loaded tungsten: A global model based on thermodynamics, kinetics and density functional theory data[END_REF][START_REF] Sun | Critical concentration for hydrogen bubble formation in metals[END_REF][START_REF] Kato | Super-saturated hydrogen effects on radiation damages in tungsten under the highflux divertor plasma irradiation[END_REF][START_REF] Ohsawa | Thermodynamics of hydrogeninduced superabundant vacancy in tungsten[END_REF]. Qualitative explanations have been proposed to explain how the kinetic barriers to form vacancies are reached during ion implantation [START_REF] Gao | Deuterium supersaturation in low-energy plasma-loaded tungsten surfaces[END_REF][START_REF] Hodille | Hydrogen super-saturated layers in h/d plasma loaded tungsten: A global model based on thermodynamics, kinetics and density functional theory data[END_REF][START_REF] Gao | Hydrogen atom-ion synergy in surface lattice modification at sub-threshold energy[END_REF] but further studies are required to establish a consistent model and so far, the creation rate used in [START_REF] Hodille | Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials[END_REF][START_REF] Ogorodnikova | Deuterium retention in tungsten in dependence of the surface conditions[END_REF][START_REF] Ogorodnikova | Ion-driven deuterium retention in tungsten[END_REF] to describe the SSL formation are ad-hoc parameters used to fit the experiments. According to the reported values, the saturation of the SSL occurs for fluence of the order of 10 24 Dm -2 which is much lower than the order of magnitude of the fluences simulated in this paper. Thus, we did not use trap creation model to simulate the SSL formation in this work.

Another aspect which is not treated here is the evolution of the SSL thanks to mono defect or defect clusters migration, clustering or annihilation with other mono-vacancies or self-interstitial atoms (SIA) created in the SSL region. The binding of H with SIA is not high [START_REF] Heinola | Hydrogen interaction with point defects in tungsten[END_REF] and the trapping of H into mono-vacancies inhibits their recombination with SIA [START_REF] Kato | Super-saturated hydrogen effects on radiation damages in tungsten under the highflux divertor plasma irradiation[END_REF][START_REF] Shutikova | The Recombination of Tungsten Intrinsic Defects near the Surface Blocked by Hydrogen and Helium Atoms[END_REF]. Thus, due to their high mobility [START_REF] Becquart | Microstructural evolution of irradiated tungsten: Ab initio parameterisation of an okmc model[END_REF][START_REF] Becquart | An object kinetic monte carlo simulation of the dynamics of helium and point defects in tungsten[END_REF], the SIA can move in the lattice, and since the SSL is very close to the surface, the SIA would most likely disappear at it. However, part of them may diffuse deeper toward the bulk, agglomerate into bigger SIA clusters which increases the binding energy of H to this type of defects [START_REF] De Backer | Hydrogen accumulation around dislocation loops and edge dislocations: from atomistic to mesoscopic scales in BCC tungsten[END_REF]. An other aspect of this SSL evolution is that empty mono-vacancies can move above 500 K [START_REF] Debelle | First temperature stage evolution of irradiation-induced defects in tungsten studied by positron annihilation spectroscopy[END_REF][START_REF] Zibrov | Vacancy cluster growth and thermal recovery in hydrogen-irradiated tungsten[END_REF]. During the steady-state of scenario I, the SSL is almost empty. The defects in the SSL could evolve and eventually form vacancy clusters containing HI [START_REF] Heinola | Modelling of the effect of elms on fuel retention at the bulk w divertor of jet[END_REF][START_REF] Zibrov | Vacancy cluster growth and thermal recovery in hydrogen-irradiated tungsten[END_REF][START_REF] Ryabtsev | Deuterium thermal desorption from vacancy clusters in tungsten[END_REF]. The detrapping energies of HI from such defects are known from DFT calculations [START_REF] Hou | Predictive model of hydrogen trapping and bubbling in nanovoids in bcc metals[END_REF]. The mobility and clustering of vacancies (or traps in general) is not yet implemented in our model but one could use either kinetic Monte Carlo [START_REF] Becquart | Microstructural evolution of irradiated tungsten: Ab initio parameterisation of an okmc model[END_REF][START_REF] Becquart | An object kinetic monte carlo simulation of the dynamics of helium and point defects in tungsten[END_REF] or cluster dynamics [START_REF] Ahlgren | Hydrogen isotope exchange in tungsten during annealing in hydrogen atmosphere[END_REF] methods to assess the distribution and evolution of the defects in the SSL to determine the trap concentration in this zone.

4.3.2.

Thermo-mechanical induced traps During scenario I and II, the temperature at the surface of the monoblock is above the ductile-to-brittle transition temperature (DBTT) of W. Because of the thermal expansion, the temperature gradient the thermal expansion mismatch between materials, residual stress could be created in the monoblock and since it is above the DBTT, plastic deformation can take place. This plastic deformation is the source of new trapping sites like dislocations but it has also been shown that strain in the 3-38% range give birth to vacancy like defects (jogs at dislocations) that are deep traps for HI [START_REF] Zibrov | Deuterium trapping by deformation-induced defects in tungsten[END_REF]. The estimation of this effect on the HI retention in the monoblock can be adressed as in [START_REF] Benannoune | Multidimensional finite-element simulations of the diffusion and trapping of hydrogen in plasma-facing components including thermal expansion[END_REF] by first computing the plastic strain induced by thermal fields, which is then used in the current modelling as a time and space dependent trap density. It is worth noting that the plastic strain estimation is strongly dependent on metallurgical features (dislocations, GB, grain sizes and orientations ...) that are triggered with temperature [START_REF] Durif | Numerical studo of the influence of tungsten recrystallization on the divertor component lifetime[END_REF].

Conclusion

In this work, the deuterium transport in W plasma facing components of the WEST and ITER divertors is simulated. Three scenarios are simulated that are relevant for plasma operations. The surface temperatures of these scenarios are 1456 K, 870 K and 435 K which are representative of temperature at the strikepoint, few cm from the strike point and far from the strike point during an attached plasma.Two types of plasma exposures are simulated: a continuous plasma exposure of 10,000 seconds and 25 cycles of 1,600 seconds with 400 seconds of steady-state plasma. Two types of simulations are run: 1D simulations in the 7.2 mm thick W armour which takes into account the super-saturated layer present in the sub-surface region and 2D simulations with the actual multi-material geometry of the components.

The 1D simulations show that there is no visible difference between the continuous and cycled plasma if one focus on the retention and migration in the bulk of the plasma facing components. This is due to moderate outgassing between discharges, which is limited by the retrapping process. The cycling introduces plasma ramp-up and ramp-down during which a dynamic retention is observed. During the ramp-down, the surface temperature decreases together with the particle flux which creates a deuterium rich layer of about 10 µm below the surface. This layer is quickly desorbed to the plasma during the following ramp-up. This dynamic retention is only visible in the hottest scenarios and is always the same independently of the cycles.

The 2D simulations are only run with the continuous plasma exposures. To speed up the simulations, the implantation is simulated with a surface concentration of mobile particles taken from the 1D simulations. From these interfaces simulations, three macroscopic quantities per monoblock are extracted: the total deuterium retention, the permeation flux to the cooling system and the desorption rate from the poloidal edge of the components. the simulations show that the retention in the super-saturated layer can reach 10 15 D per monoblock (for the coldest exposures) but for long plasma exposure, the retention is dominated by the bulk. However, the presence of such high quantity of D near the surface may affect the edge plasma, especially during transient. Then, as already mentioned in previous studies [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF][START_REF] Hodille | Estimation of the tritium retention in ITER tungsten divertor target using macroscopic rate equations simulations[END_REF][START_REF] Sang | Modeling of fuel retention in the upper tungsten divertor of east from attached to detached divertor plasma[END_REF] despite being the scenario with the highest incident flux of particles, the hottest scenario has the lowest retention. It is also the only scenario leading to the permeation of D atoms after 10,000 seconds of exposure. Finally, it is calculated that the desorption from the toroidal edges of the components stays a negligible source of recycling depending on the morphology of the heat and particle flux on the divertor.

3. 1

 1 .1. D depth profile Figure3shows the evolution of the D depth profile in the 7.2 mm W armour after 5 -25 cycles for scenario I (a), II (b) and III (c). The depth profile just before the ramp-down phase (end of the steady-state plasma) is also shown (black dashed lines). The D depth profiles are similar in scenarios I and II after the end of the cycles: a high concentration in the SSL area (5-10 nm) with about 10 27 D/m -3 , a second step up to 10 µm with about 10 25 D/m -3 accompanied by a tail up to 1 mm and the bulk part. The D depth profiles for scenario III only exhibit two features: the high concentration in the SSL and the step of constant concentration up to the migration depth, which after 25 cycles is around 400 µm.

10 - 10 Figure 3 :

 103 Figure 3: Results of the MHIMS simulations: Concentration of D in the W armour after 5 -25 cycles of 1,600 seconds with the exposure conditions of scenario I (a), scenario II (b) and scenario III (c). The concentration of D before the ramp down of cycle 25 is also shown in black dash line.

- 2 )Figure 4 :

 24 Figure4: Evolution of the D retention in the 7.2 mm W armour after 10,000 seconds of steady-state continuous plasma (dashed lines) or after 25 cycles (solid line). The retention is shown as a function of the implanted fluence to be able to compare continuous and cycles exposures.

Figure 5 :

 5 Figure 5: Evolution of the D retention with time during the 10 th cycle of scenario II (T SSP surf = 870 K. The beginning/end of the plasma ramp-up/down are marked with vertical dashed line and the retention released during the ramp-up (RU rel) and retained during the ramp down (RD ret) are marked with arrow.

  ret Sc. II, RD ret Sc. III, RD ret Sc. I, RU rel Sc. II, RU rel Sc. III, RU rel

Figure 6 :

 6 Figure 6: Quantity of D retained during the ramp down (RD) and released during the ramp-up (RU) phases as a function of the number of cycles N cycle for the different scenarios.

Figure 7 :

 7 Figure 7: Results of the FESTIM simulations for scenario I (a) and II (b): deuterium retention field c D (m -3 ), in the semi-monoblock after 10,000 seconds of continuous. The colour bar have different scales .

3. 3 .(

 3 Macroscopic quantities 3.3.1. D inventory The total D inventory in one monoblock exposed to scenario I, II or III can be calculated as: Dm -1 ) is the retention computed by FESTIM in the 2D geometry which accounts for the intrinsic traps and the mobile particles. R SSL,MHIMS D (Dm -2 ) is the retention computed by MHIMS in 1D in the SSL only, as the mobile particles and the HI trapped in the intrinsic traps are already taken into account in R FESTIM D

Figure 8 :

 8 Figure 8: Results of the FESTIM simulation for the scenario III: deuterium retention field c D (m -3 ) after 10,000 second of continuous plasma exposure (on Γ top ).The thickness of the simulated material is 500 µm. The width is the same as for the other scenarios (14 mm) but for sake of clarity, the figure is zoomed in the corner from which the D is allow to desorb (Γ lateral ). The temperature is 435 K.

Figure 9 :

 9 Figure 9: D retention per monoblock I D for the three simulated scenario during a continuous 10,000 s plasma exposure. The contribution of the SSL retention is also shown with the dashed line.

Figure 10 :

 10 Figure 10: D release to the cooling pipes(a) and from the toroidal edges (b) extracted from the 2D-FESTIM simulations for scenario I, II and III. For scenario II and III, the release to the cooling pipes (permeation flux) is 0 and both curves overlap each other.

Table 4 :

 4 Output value of c MHIMS m (X imp ) given as input to the 2D FESTIM simulations. This is the local concentration of mobile particles in the implantation zone.

	scenario	c MHIMS m	(X imp ) (m -3 )
	I II III	2.01 × 10 22 3.20 × 10 22 9.74 × 10 22

Acknowledgements

Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the EURATOM research and training program 2014-2018 and 2019-2020 under grand agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The project leading to this publication has received funding from the French National Research Agency (Grant No. ANR-18-CEO5-0012).