Parametric study of hydrogenic inventory in the ITER divertor based on machine learning

A parametric study is performed with the 2D FESTIM code for the ITER monoblock geometry. The influence of the monoblock surface temperature, the incident ion energy and particle flux on the monoblock hydrogen inventory is investigated. The simulated data is analysed with a Gaussian regression process and an inventory map as a function of ion energy and incident flux is given. Using this inventory map, the hydrogen inventory in the divertor is easily derived for any type of scenario. Here, the case of a detached ITER scenario with inputs from the SOLPS code is presented. For this scenario, the hydrogen inventory per monoblock is highly dependent of surface temperature and ranges from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{18}$$\end{document}1018 to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$6 \times 10^{19}$$\end{document}6×1019 H after a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{7}$$\end{document}107 s exposure. The inventory evolves as a power law of time and is lower at strike points where the surface temperature is high. Hydrogen inventory in the whole divertor after a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{7}$$\end{document}107 s exposure is estimated at approximately 8 g.

Scientific Reports | (2020) 10:17798 | https://doi.org/10.1038/s41598-020-74844-w www.nature.com/scientificreports/ Methodology Model description. As described in previous studies 5,[7][8][9] , the macroscopic rate equations model used in this work splits hydrogen isotopes into two populations: the mobile particles and the trapped ones. The temporal evolution of mobile particles c m and trapped particles c t,i in the i-th trap are described in Eqs. (1) and (2) respectively.
In Eq. (1), D(T) = D 0 · exp − E diff /(k B · T) is the diffusion coefficient in m 2 s −1 , T the temperature in K and k B = 8.617 × 10 −5 eV K −1 the Boltzmann constant, S is the volumetric source term of mobile particles in m −3 s −1 (which can take into account plasma implantation), k(T) = k 0 exp − E k /(k B · T) and p(T) = p 0 exp − E p /(k B · T) are the trapping and detrapping rates expressed in m 3 s −1 and s −1 respectively. n i is the trap density in m −3 . Concentration conservation is assumed at interfaces for the sake of simplicity and computation time. It has already been shown that retention in monoblocks is dominated by retention in tungsten 5 . Therefore, conservation of chemical potential will not have much influence on the results. A more complete description of this model is given in 5 .
The temperature temporal evolution is governed by the heat equation described as follow: where ρ is the density of the material in kg m −3 , C p its specific heat capacity expressed in J kg −1 K −1 and the thermal conductivity expressed in W K −1 . Equations (1), (2) and (3) are then solved in FESTIM using the finite element method implemented in the FEniCS project 10 . FESTIM is implemented in Python and provides a user-friendly interface for performing multiphysics, multidimensional and multi-material simulations 5 . All plots in this work were generated with Matplotlib 11 . Simulation description. The first step of the work was to simulate the hydrogenic transport and trapping in a tungsten monoblock as a function of the loading conditions. Moreover, a parametric study will be carried out in order to simulate the whole range of the implantation conditions encountered in the ITER divertor.
Geometry. The geometry used in this work is that of a non-shaped ITER monoblock (see Fig. 1). The monoblocks use tungsten armour and a 1.5 mm-thick CuCrZr pipe as heat sink. The pipe is jointed to the tungsten. A 1 mm-thick Cu interlayer is used in order to handle stress resulting from differential thermal expansion 12 . The surface Ŵ top is facing the plasma and Ŵ coolant is cooled by water.
Material properties. The material properties used in these simulations are described in Table 1 and their temperature dependence is shown in Fig. 2. The trap parameters are described in Table 2. Influence of mechani- cal fields such as thermal expansion on trap creation 7 was not taken into account in this work. Hodille et al. described an extrinsic trap in tungsten created by ion implantation 9 . This trap is assumed to have only a small influence on the macroscopic behaviour of the monoblock and is therefore not taken into account in this work for the sake of simplicity.
Boundary conditions. Mobile particles concentration c m is imposed on Ŵ top which allows to simulate particle implantation without having to include a volumetric source term applied on the first few nanometres. This approximation allows to have a broader mesh and therefore saves computation time without affecting the macroscopic behaviour. Molecular recombination is assumed on Ŵ coolant . Even though it could be assumed on the gaps between monoblocks, it can be shown that its influence on the macroscopic behaviour remains low. Desorption from the other surfaces is therefore assumed to be zero for simplification purposes. Uniform heat loads ϕ H are applied on the surface Ŵ top with a Neumman boundary condition or temperature is constrained on Ŵ top with a Dirichlet boundary condition and a convective exchange condition is set on surface Ŵ coolant . All the other surfaces are assumed thermally insulated. The set of boundary conditions can finally be described as follow: Table 1. Materials properties used in the simulations. Thermal properties are fitted from ANSYS [13][14][15] .

Results
Thermal behaviour. Steady-state heat transfer simulations were performed with FESTIM with varying heat loads ϕ H . With ϕ H = 1 MW m −2 , the surface temperature of the monoblock was found to be around 400 K (see Fig. 3a) whereas with ϕ H = 10 MW m −2 the surface was around 1400 K (see Fig. 3b).
In order to simplify the analytical relations used in "Influence of incident particle flux and ion energy on hydrogen inventory" section, only the mean surface temperature was considered in the following sections. T surface therefore increases linearly with the heat load and can be modelled by Eq. (8) (see Fig. 3c).
This was found to be in very good agreement with experimental measurements performed in 18 .   The assumption of a constant surface temperature had low influence on the results compared to a non-homogeneous surface temperature that could be obtained with a heat flux condition since surface temperature gradient was low compared to the one between the top surface and the cooling surface. For surface temperatures below 500 K, 1D simulations were performed for the penetration depth of hydrogen remained very low (a few microns) and 1D approximation was sufficient 19 . For temperatures above 500 K for which edge effects become dominant, 2D simulations have been performed. After 10 7 s a high retention zone appeared far from the exposed surface Ŵ top (see Fig. 4). As described in 5 , this is due to thermal effects. As seen in Fig. 3a and b, the temperature was found to decrease in regions close to the cooling pipe Ŵ coolant leading to an increase in trap occupancy, creating this high retention zone. This was however not true for monoblocks where T surface ≈ T coolant since the temperature gradient in the domain is very low. Instead, trap occupancy was close to one and the retention was high in the whole region where hydrogen had penetrated and not only far from the top surface.
Hydrogen inventory in monoblocks as a function of T surface and c surface is shown in Fig. 5. In order to obtain this continuous field, more than 600 simulations randomly distributed on the parameter plane were run and analysed using a Gaussian process machine learning algorithm 20 as in 21 based on the python package inferencetools 22 . In Fig. 5, the inventory obtained by the Gaussian regression process is given for a constant value of c surf = 2 × 10 21 m −3 (top inset) and a constant temperature T = 850 K (left inset). The Gaussian regression process was particularly appropriate as it calculates a confidence interval based on the standard deviation σ . As expected, the lower the density of simulation points, the higher was the value of σ (for example around 850 K on the top inset of Fig. 5). However, despite the lack of simulation in this region, the value of σ was still acceptable (only a few percents of the inventory) ensuring the quality of the resulting interpolation.
As expected, inventory was found to globally increase with c surface . For T surface > 550 K , the inventory tended to decrease with surface temperature. However, for T surface < 550 K , inventory increased with surface temperature. This phenomenon is due to a trade-off between an increase of the detrapping rate and an increase of the diffusion coefficient making the hydrogen particles penetrate deeper into the bulk. Above 550 K , detrapping becomes dominant and inventory decreases. This mapping of inventory as a function of T surface and c surface provides an easy way of estimating the inventory in monoblocks for several exposure conditions without having to run many simulations. Indeed, to estimate the inventory with different exposure conditions, one only needs to associate these conditions (ϕ inc , E) to a couple (c surf , T surf ).
Influence of incident particle flux and ion energy on hydrogen inventory. Incident particle flux ϕ inc and ion energy E have an impact on the amount of mobile particles implanted in the material but also on the heat load and therefore on the surface temperature of the monoblock. www.nature.com/scientificreports/ Assuming a source term with a narrow Gaussian distribution and a non-instantaneous recombination (characterised by a recombination coefficient K), the concentration c max at the near surface is approximated by: where ϕ imp = (1 − r) · ϕ inc is the implanted particle flux, r is the particle reflection coefficient, K is the recombination coefficient, and R p is the mean implantation depth in m. Details can be found in "appendix" as Supplementary Material. Many different values of the recombination coefficient K for tungsten can be found in literature. For instance the widely used Anderl coefficient describes an endothermic recombination 23 whereas Ogorodnikova showed an exothermic recombination coefficient could be used to reproduce a set of experiments 24 .
Facing the difficulty of an accurate choice for K and following the recommendation of Causey et al. 25 , an instantaneous recombination will therefore be assumed (i.e. K → +∞ ). It is also worth noting that experiments by Bisson et al. 26 support the fact that recombination is not the rate limiting step during the hydrogen release from polycrystalline tungsten after ion implantation.
In the following, the concentration on Ŵ top was set to c surface = c max for the kinetics involved are really fast (see appendix of 27 ) and R p is small compared to the monoblock dimensions.
The heat load was assumed to evolve as a function of the incident particle flux ϕ inc and E as follow: with e = 1.6 10 −19 C . This relation was obtained by fitting SOLPS data 28,29 . The factor 2.2 was applied to take into account other heat sources such as radiative flux. Moreover, the ion energy E has an influence on r and implantation range R p and it was possible to model the evolution of these parameters with SRIM 30 calculations as follow:  Fig. 6a. From the thermal behaviour given by Eq. (8), the surface temperature T surface can be computed (see Fig. 6b). Finally, c max was obtained from Eqs. (9) and (12) (see Fig. 6c). One must be aware that above 1500 K, W recrystallisation can occur and H transport will strongly be affected. The hypothesis made above as well as material properties may then not be valid. Because of the trade-off between the amount of implanted particles and the resulting heat flux, the maximum value of c max was found to be 2 × 10 22 m −3 around (ϕ inc , E) = (8 × 10 22 m −2 s −1 , 20 eV) . Considering the previously calculated response of the monoblock to c surface and T surface (see Fig. 5), the inventory as a function of ϕ inc and E was computed (see Fig. 6d). The inventory values have not been calculated for surface temperatures above 1200 K. Again a trade-off was found between implanted particle flux and surface temperature. Indeed, the maximum inventory was not found at regions where the incident flux is maximum but rather at regions where c surface is maximum and T surface is minimum as seen in previous 1D studies 8 . Application to tokamak exposure conditions. Each white circle in Fig. 6 corresponds to a point along a poloidal section of the ITER divertor for which implanted particle flux and ion energy were calculated with SOLPS 6 for a partially detached plasma scenario. This scenario corresponds to a Q = 10 discharge with a neutral pressure of 8.6 Pa 31 .  Figure 6. ϕ H , T surface , c max and inventory per monoblock as a function of ϕ inc and E. Inventory has not been calculated for surface temperature above 1200 K (greyed region). White circles correspond to points on ITER divertor using the divertor plasma parameters from SOLPS 6 calculations (see "Application to tokamak exposure conditions" section).
Scientific Reports | (2020) 10:17798 | https://doi.org/10.1038/s41598-020-74844-w www.nature.com/scientificreports/ As expected the highest surface temperatures and heat loads were located on strike points and most of the zones on the divertor were found to stay at coolant temperature (see Fig. 7). The maximum hydrogen content is approximately 6 × 10 19 H per monoblock after a 10 7 s exposure. As explained in the previous section, the maximum inventory is not necessarily in the region where the flux is maximum as it induces a higher temperature which will tend to increase detrapping: strike points are not where hydrogen is trapped the most. Instead, the maximum inventory is reached about 5 cm away from the strike points where the temperature and the fluxes are high enough to guarantee a strong source of mobile particle but the temperature is not high enough to trigger detrapping.
For all points on the divertor, the inventory evolved as a · t b as shown in Fig. 8 for particular points on the inner vertical target ( x = 0.03 m is close to the strike point). The coefficient b is maximum on strike points reaching 0.75 (see Fig. 9). In other regions, b is closer to 0.5. This result can be explained by the non-homogeneous temperature field in monoblocks with high heat loads. For monoblocks with a high surface temperature, as hydrogen penetrates deeper into the bulk, the bulk temperature decreases (see Fig. 3b) leading to an increase of the trap occupancy 8 . The exponent b is therefore higher than 0.5. For monoblock where T surface ≈ T coolant on the other hand, the temperature is homogeneous in the whole domain and b = 0.5 . This corresponds to a diffusion-limited behaviour.
The temperature is close to T coolant and the trap occupancy is therefore close to one in the whole domain which is not the case for regions near strike points where temperature fields are non-uniform.
One can obtain the inventory in the whole divertor by integrating the results obtained in Fig. 7 over the tokamak as follow: with N cassettes = 54 the number of cassettes, N PFU−IVT = 16 and N PFU−OVT = 22 the number of plasma facing units per cassette on the inner and outer targets respectively, inv IVT and inv OVT the hydrogen inventory profile along the inner and outer targets respectively and x the distance along the targets.
After a 10 7 s exposure, hydrogen inventory is estimated at approximately 8 g (see Fig. 10) which is relatively low considering the ITER in-vessel limit and the elapsed time. De Temmerman et al. showed that retention in ITER can reach 0.3 g per 400 s discharge when taking into account Be deposits. First, a steady state exposure was considered for simplification purposes. This result is however conservative. As seen in 5,8 , cycling effects could have an influence in regions where T surface varies a lot, for example within 10 cm on both sides of the strike points. Though, since a large majority of monoblocks were found to stay at room temperature, even during operations (see Figs. 6b and 7) the thermal effect should remain low and discrepancies would rather be due to particle flux evolution along the target.
Shaping of monoblocks (e.g. chamfers) was not taken into account in this work for simplification purposes. Such shaping can have an influence on the incident particle and heat loads on the plasma facing surface of the monoblocks.
This study presents the hydrogen trapping in W monoblocks. It shows that the latter remains low but, as already pointed out by JET studies, the trapping on Be co-deposited layers is expected to be the main mechanism for tritium retention in ITER 32,33 . Such layers could be found in the cold regions of the divertor but as soon  www.nature.com/scientificreports/ as the strike points hit these layers, they should be sputtered away (as sputtering of Be is possible even at low energy 32,34 ). The retention where the deposited layers are not present (either sputtered or not formed anyway) would then be given by the model presented here. The molecular recombination coefficient at the surface of the cooling pipe was taken from 17 and was measured in vacuum. One could argue that recombination in presence of water will be facilitated. It can however be shown that this parameter has very low influence on the inventory since it was dominated by retention in tungsten. This parameter will however have an influence on the permeation flux and should be studied in future work.
Similarly, the influence on molecular recombination on the sides of the monoblock was found to have a low impact on the results. By assuming an instantaneous recombination coefficient, the relative error on the monoblock inventory was found to be significant only in hot regions (i.e. within 10 cm on both sides of the strike points). The influence on the total divertor inventory is therefore low (less than 5% after a 10 7 s exposure) since it is dominated by regions where T surface ≈ T coolant .
It should be noted that specific scenarii like edge localised modes (ELMs) were also not taken into account in this work since their time scale is very short. MRE simulations by Hu and Hassanein 35 suggest that a 400 s discharge with 1 Hz or 10 Hz ELMs significantly reduces (77%) the inventory in W materials. However, the modelling of the ELM is simulated by increasing the temperature for a very short time without changing the incident flux of particles that can also be much higher thus balancing the fuel retention reduction. Another study by Schmid et al. 36 also simulated the effect of 1 Hz ELMs on fuel retention in W. The outcome is that 6 s of 1 Hz-ELMs does not affect significantly the fuel retention, though the temperature excursion in those simulations are smaller than for the one of Hu and Hassanein. Thus, the effect of ELMs, especially the balance between increase of heat flux, incident energy and particle flux, could either favour or disfavour trapping, diffusion and migration and therefore the overall retention.
In this study the model to link the concentration of mobile particles at the surface (implantation zone) with the exposure condition considers that the particles are implanted in the bulk and that the recombination coefficient is very high since many uncertainties concerning the recombination coefficient are yet to be lifted. However, if an exothermic process is considered as in 24 , this should have low influence since recombination is very quick at a temperature close to that of the coolant.
On the other hand, experimental results 37 suggest that for ion energy below 5 eV/H, typical of detached plasma as the one treated in the previous section, the surface process can be important and limits the uptake of hydrogen, i.e. the adsorption on the surface and the further absorption from surface to bulk could be the limiting process for the growth of c surface during such exposure. The evolution of c surface to the exposure condition for that range of energy would therefore be different and therefore the inventory. The advantage of the presented method is that taking into account such process is relatively easy as no expensive simulations are needed. One would only need to modify the model giving c surface as a function of (E inc , ϕ inc ) to include the different surface processes. To this end, one can use kinetic surface models [38][39][40][41] .
Trap properties have a great impact on the inventory. In this study, a homogeneous trap distribution is assumed for simplification purposes. A more thorough study could investigate the influence on trap distribution, energy and density. Trap properties might also vary along the divertor based on exposure conditions. Moreover the impact of neutrons must be assessed as neutron-induced traps have a high detrapping energy.
Finally, helium implantation in the materials and bubble formation could modify the hydrogen transport in monoblocks.

Conclusion
ITER-like monoblocks have been studied using a novel method in order to estimate the hydrogen content as a function of exposure conditions such as implanted particle flux, ion energy, heat load and monoblock surface temperature. Several hundred data points have been simulated with FESTIM and analysed to estimate the hydrogen inventory in monoblocks for any input conditions using Gaussian regression, a machine learning algorithm www.nature.com/scientificreports/ which calculates the confidence interval for each point, Thanks to this relation, one can easily estimate hydrogen content in the whole divertor without having to run all the simulations. An application has been made based on the output from a SOLPS calculation of exposure conditions distribution on the ITER divertor and shows that for these conditions the inventory could reach 10 20 H per monoblock near strike points after a 10 7 s exposure. The total hydrogen content in ITER divertor is estimated to be 8 g which is well below the inner-vessel safety limit of 1 kg. Future work will include applying this technique to calculations performed on the WEST tokamak and studies of the impact of trap parameters (especially neutron-induced traps) on the inventory results as well as the impact of more complex geometries.