Impact of Soret effect on hydrogen and helium retention in PFC tungsten under ELM-like conditions

In our previous work, we have demonstrated using nonequilibrium molecular-dynamics simulations that the fluxes of helium and self-interstitial atoms in the presence of a thermal gradient in tungsten are directed opposite to the heat flux, indicating that species transport is governed by a Soret effect, namely, thermal-gradient-driven diffusion, characterized by a negative heat of transport that drives species transport uphill, i.e. from the cooler to the hot regions of the tungsten sample. In this work, the findings of our thermal and species transport analysis have been implemented in our cluster-dynamics code, Xolotl, which has been used to compute temperature and species profiles over spatiotemporal scales representative of plasma-facing component (PFC) tungsten under typical reactor operating conditions, including extreme heat loads at the plasma-facing surface characteristic of plasma instabilities that induce edge localized modes (ELMs). We demonstrate that the steady-state species profiles, when properly accounting for the Soret effect, vary significantly from those where temperature-gradient-driven transport is not accounted for and discuss the implications of such a Soret effect on the response to plasma exposure of plasma-facing tungsten. Although our cluster-dynamics simulations do not yet include self-clustering of helium or hydrogen blister formation, our simulation results show that the Soret effect substantially reduces helium and hydrogenic species retention inside PFC tungsten.


Introduction: transient heat load under edge localized mode (ELM)-like conditions
Under burning plasma conditions, plasma-facing components (PFCs) in the International Thermonuclear Experimental Reactor (ITER) and future fusion reactors with magnetic confinement are expected to experience harsh conditions. For example, the divertor wall of ITER will experience a high particle flux of ∼10 20 -10 24 m 2 s −1 and a time-averaged thermal load up to 20 MW m −2 under steady-state operating conditions. In addition, during transients like type I ELMs in the H-mode plasma confinement regime, a very high, repetitive thermal load up to ∼1 GW m −2 are periodically released from the plasma over a very short time scale, typically 100-500 µs. Under reactor-relevant conditions, individual particle fluxes will cause significant damage to the tungsten divertor target. For example, self-clustering of implanted helium will damage the near-surface region of PFC tungsten (main candidate as plasma-facing material) and, in turn, drive the tendril-like growth of surface features; hydrogen and hydrogen isotopes will lead to blister formation and also result in embrittlement of the wall. In addition to these impurity species, 14 MeV neutrons generated in the fusion reaction will penetrate the PFC material creating intrinsic defects, such as vacancies, self-interstitial atoms, and clusters of such point defects. Furthermore, defect clusters not only compromise the mechanical strength of the PFC materials but also affect the thermal conductivity of PFC tungsten as has been shown by both experimental observations [1,2] and molecular-dynamics simulations [3]. Unmitigated ELMs in ITER coupled with low thermal conductivity in the helium-implanted near-surface PFC tungsten will lead to surface melting or cracking [4], which would compromise the structural integrity of the divertor target tungsten monoblocks.
In our previous study [5], we have demonstrated using non-equilibrium molecular-dynamics (NEMD) simulations that intrinsic point defects and impurities in tungsten show thermal-gradient-driven diffusion, also known as Soret diffusion, characterized by a negative heat of transport that drives species transport uphill, i.e. from the cooler to the hot regions of the tungsten sample. In this work, we study the species transport, focusing, in particular, on helium, deuterium, and tritium, under an ELM-like heat load to understand how the Soret effect influences hydrogen and helium retention, which have a significant impact on the mechanical strength of the divertor wall, and tritium inventory, which is critical from a safety perspective [6,7]. Implications of the Soret effect on the dynamical response of PFCs in fusion reactors has been previously studied [8,9] but not under extreme ELM-like heat loads, where the Soret effect is expected to be significant as a result of the strong temperature gradients anticipated to be developed. Our current work systematically addresses this question.
The rest of the article is structured as follows. The models for heat and species transport, the relevant model parameters, and the simulation method used in our study are described in section 2. The main results of our simulations are presented in section 3, culminating with the impact of the Soret effect on species retention. A succinct summary and a discussion of our results and their implications for PFC tungsten response to plasma exposure are given in section 4.

Simplified tungsten monoblock geometry
Tungsten is the candidate material for the inner and outer divertor walls of ITER due to its mechanical strength, high melting point, low anticipated tritium retention [10,11], and many other thermophysical properties which are favorable for tolerating the harsh environment expected in the divertor region [12]. The divertor will be fabricated using tungsten monoblocks [13][14][15] that have dimensions of 28 mm × 26-28 mm × 12 mm with a CuCrZr cooling tube of inner diameter of 12 mm and thickness of 1.5 mm located at the monoblock center, and mounted with a Cu interlayer of thickness 1 mm in between, as illustrated in figure 1(a). The minimum tungsten armor thickness in the design (corresponding to the shortest distance from the monoblock top surface to the copper interlayer) is 6 mm.
The purpose of this study is to obtain a fundamental understanding of how the high thermal gradients generated by ELMlike heat loads affect the transport of atomic species, such as helium and hydrogenic atoms and intrinsic point defects, and, in turn, influence the retention of helium and hydrogenic species in the PFC tungsten. Towards that goal, we have simplified the monoblock geometry to consist of three slabs with thicknesses of 6 mm, 1 mm, and 1.5 mm, made of tungsten, copper, and CuCrZr, respectively, as shown in figure 1(b). We will later show that the 1D temperature profile along the depth of the slab, which is obtained by solving the heat balance equation under the appropriate initial and boundary conditions, for the simplified slab-like geometry, compares very well with the predicted profiles using the full 2D monoblock geometry.

The ELM-heat-load model
The energy balance equations for 1D heat transport in each of the slabs in the simplified monoblock geometry are expressed by where T(z, t) is the evolving temperature field and ρ s , C v,s , and κ s are the density, specific heat, and thermal conductivity of the constituent materials in each slab, i.e. tungsten, copper, and CuCrZr alloy. In the simplified geometry, the three slabs are connected in series and continuity of temperature, as well as continuity of heat flux at the tungsten-copper and copper-CuCrZr interfaces are assumed; in other words, the two interfaces are assumed to be in thermal equilibrium (equal temperature at each side of the interface) and there is no heat production or loss at either interface, i.e. the heat flux that arrives at each interface side is equal to the heat flux conducted at the other side of each interface. Additionally, the  temperature dependence of the specific heat and thermal conductivity of the slab materials [16] is accounted for. For the boundary condition at the plasma-facing side of the tungsten slab, we have used an empirical model for the ELM temporal heat load profile to express the total perpendicular heat load under ELM-like conditions, q ⊥,total (t) [16]. We have neglected the spatial variation of the heat load over the plasma-facing surface and assumed that the spatial peak of the unsteady ELM heat load q ⊥,ELM (t) and the inter-ELM average heat flux q ⊥,interELM are occurring at that same location with these two loads together constituting the total heat load as with τ = 0.8τ IR , where τ IR is the time constant of the rise phase of the ELM. To simplify our calculations, we have introduced a new quantity q eqv , representing an equivalent heat load, which is defined as where τ ELM is the ELM period (i.e. the inverse of the ELM frequency). A representative heat-load profile under ELM-like conditions is shown in figure 2(a). For the boundary condition at the coolant side of the monoblock model, we have used a set of empirical correlations, equation (4), to express the heat removal q loss by forced convection where the coolant liquid is water at 100 • C under 30 bar pressure and flowing with a velocity of 12 m s −1 through a swirl tape with thickness of 2 mm and twist ratio of 2, as specified in [16,19]. Equation (4) is the best-fit of the heat-loss data as a function of coolant-side temperature under the aforementioned coolant conditions, obtained through [17].
The raw data and the optimally fitted empirical relations are plotted in figure 2(b). Under these coolant conditions, the heat loss is below the critical heat flux which is approximately 61 MW m −2 according to the Tong-75 correlation [18].

Modeling predictions for the monoblock geometry and comparison with those employing the composite slab and tungsten slab approximations
The evolution of the temperature under heat load, simulated using a 2D monoblock geometry similar to that in figure 1(a), has been reported in the literature [16,20]. While the absolute value of the maximum temperature reported in [16] under a heat load identical with that used in this work is ∼400 K higher than the maximum temperature predicted using the simplified composite slab-like geometry, the predicted difference between the maximum and the minimum surface temperature is the same in both cases. However, under the different heat-load condition studied in [20], the 2D temperature profile shows a significant temperature variation in the lateral direction at a fixed depth from the plasma-facing surface and the predicted maximum temperature is located at the corner of the monoblock. Therefore, it would be reasonable to assume that the predicted 1D profile in our simplified slab geometry is much closer to an average temperature profile, with the averaging performed over the lateral dimension.
Furthermore, it is evident from figure 2(c) that the temporal variation in the temperature profile (as measured by the difference between maximum and minimum temperature) for the composite slab-like geometry is limited to within 1 mm (b) Dependence of heat loss due to forced convection on the coolant-side wall temperature as predicted by the Nukiyama curve [17,18] (blue solid line) together with a best-fitted correlation obtained from the data (red dashed line). (c) 1D temperature profile across the slab when the plasma-facing surface temperature is at its maximum and at its minimum under the ELM-like heat-load cycle and for two different types of slabs: a composite slab with a thickness of 8.5 mm (as shown in figure 1(b)) and a pristine tungsten slab with a thickness of 6 mm, and under qeqv heat load for a pristine tungsten slab with a thickness of 6 mm. from the plasma-facing surface and any variation beyond that depth is insignificant. Therefore, the simulation of the species concentration evolution can be further simplified by restricting the simulation domain to a tungsten slab with a thickness of 6 mm and a modified thermal flux at the coolant side. In figure 2(c), the resulting temperature profile is plotted in orange, which is practically identical to the profile predicted when employing the composite-slab approximation and, hence, justifies this additional simplification. In this study, we have not accounted for self-clustering of helium that results in formation of over-pressurized helium bubbles and, therefore, decrease in the thermal conductivity of the tungsten. However, in the simulations, we have included the decrease in the thermal conductivity of tungsten due to helium implantation in a thin near-surface region with a thickness of 190 nm, a depth consistent with the choice of the energy of the implanted particles, and termed this special condition as 'damaged tungsten' in order to distinguish it from the pristine tungsten results without considering this (thermal conductivity reduction) effect. Finally, under a fixed heat load defined by equation (3), equivalent to the time-periodic ELM-like heat load, the temperature profile matches closely with the 6 mm pristine tungsten slab, as well as the composite slab except for the first 1 mm from the plasma-facing surface. Despite this difference in the temperature profile, a marginal difference between the total species content predicted under equivalent and ELM-like heat load along with the significant computational advantage of using a fixed equivalent heat load makes this equivalent heat load a more practical approach to simulate ELM-like conditions, which is discussed later.

Species transport-Soret effect
The species balance equations for individual species along with the appropriate boundary and initial conditions are used to simulate the evolution of the distribution of helium, deuterium, and tritium in tungsten. The species concentration evolution due to a combined effect of Fickian diffusion and thermo-migration (i.e. Soret diffusion) is governed by [5] where the diffusion coefficient D A of species A in tungsten is expressed through the Arrhenius equation with D 0 A and E A being the pre-exponential factor and activation energy barrier for diffusion, respectively, and the source term Γ A is the volumetric species implantation profile which is strongly dependent on the irradiation conditions. The boundary conditions for the species concentration fields are free surfaces (zero species concentration) on both sides of the material, as it is sufficient for typical laboratory implantation conditions [21].

Model parameters
Specific heat capacities and thermal conductivities for the different materials are taken from table B.1 of [16]. The parameters used in the equations of section 2 are listed in table 1. An incident flux of 10 21 m −2 s −1 is used for all three specieshelium, deuterium, and tritium-with an implantation profile for each species obtained for pristine tungsten through binary collision approximation [22] simulations corresponding to ITER burning plasma operation at peak heat flux in the divertor as described in [23]. In this study, the incident species flux is taken to be constant, as the main focus of this work is to examine the impact of the Soret effect. The time dependence of the incident species fluxes will be accounted for in subsequent studies, with further complexity added to our simulations. We have simulated two different situations: a pristine tungsten slab and a damaged (helium implanted) tungsten slab, both with a thickness of 6 mm. The damaged tungsten slab is modeled as consisting of a thin layer of thickness of 190 nm near the surface with a lower thermal conductivity than that of pristine tungsten (equal to 20% of that of pristine tungsten), with the remainder of the material in the slab considered as pristine.
The Soret coefficient (heat of transport) for hydrogen is calculated following the same method used in [5] to calculate the Soret coefficient for helium. In this study, we have assumed that the Soret coefficients for deuterium (Q * D ) and tritium (Q * T ) have the same values, equal to that of hydrogen; however, the diffusivities of the various hydrogen isotopes are taken to be mass dependent as listed in table 1.

Simulation methods and their implementation in Xolotl
Our cluster-dynamics code, Xolotl [24], solves well-posed boundary-value problems for the concentration fields of the various species in an irradiated material, here tungsten, governed by the corresponding spatially dependent driftdiffusion-reaction equations subject to the appropriate initial and boundary conditions. We use a cell-centered finitedifference implementation in which the grid along the depth direction has a variable grid size, with a finer grid resolution near the surface, in order to better describe the species concentration profile of the incoming impurity flux. The PETSc solver [25] is used and allows for message passing interface (MPI) parallelism by splitting the grid (discretized domain) across processors. To aid with convergence, the grids used to solve the species transport equation (5) and heat equation (1) do not need to be identical but use the same number of grid points; in this study, the grid used for the representation of the temperature field is coarser near the surface than that used for representing the species concentration fields. The temperature values are linearly interpolated between the two grids to be used in solving the species transport equations.
In the case of equivalent heat load, the steady-state temperature profile (along depth z) is first obtained by solving the heat equation. To compute the temperature profile, the transient heat conduction equation is solved until a steady state is achieved. A new simulation is then started solving for the species concentration profiles using the steady-state temperature profile that is previously obtained. We are able to separate the two equations (heat equation from species balances) because the heat equation can be decoupled from the species equations as explained in [5]. Therefore, we have independently solved the heat balance equation to obtain the evolving temperature profile, T(z, t), (see figure 3) which then is used, together with its gradient, in the species balance equations to solve for the species concentration profiles. Figure 3 shows the surface temperature evolution for three cases: pristine tungsten under an ELM-like heat load, damaged tungsten under the same ELM heat load, and pristine tungsten under an equivalent constant heat load as defined in equation (3). The temperature results without taking the Soret effect into account are not shown as they are identical to those with the Soret effect accounted for. In all cases, a steady periodic temperature state is reached within 2.1 s of exposure to the heat load (99% of the mean value of the maximum and minimum temperature levels in the oscillatory surface temperatures on the way to achieving the steady periodic state is reached after 2.1 s); after this initial transient, the amplitude and period of the surface temperature (time periodic function) remain constant. The inset shows a magnification of the surface temperature evolution at the steady periodic  [16], except for the amplitude of the time periodic state reached due to the fact that we use a simplified 1D model (of the 2D monoblock model) as discussed in section 2. The effect of the 190 nmthick layer of damaged W near the surface is not visible on the scale of the figure but leads to a slight increase of the maximum surface temperature from 1647 K, in pristine tungsten, to 1650 K. In the case where an equivalent constant heat load is used, there is no cyclic variation of temperature with time, as expected, and the surface temperature reaches a steady-state value of 1316 K (black dot-dashed line in figure 3). This temperature steady state is reached within 2.2 s of exposure to the heat load, similar to the ELM-like simulations.

Species content evolution
3.2.1. Helium. The helium content integrated over the full depth of the tungsten slab is shown as a function of time in figure 4(a). Considering the small effect of the 190 nm nearsurface layer of damaged tungsten on the temperature behavior, it follows that the helium content is only marginally affected by the presence of this thin near-surface layer as well (see values in table 2). As a result of this minimal difference, the concentration profiles of helium in damaged and pristine tungsten slabs under identical heat-load conditions are practically indistinguishable. When the Soret effect is not accounted for, the overall helium content in the tungsten slab is ∼3.5 times higher than when the Soret effect is accounted for. This substantial impact of the Soret effect on the helium content can be explained by the additional transport term (Soret diffusion) that moves helium atoms from the colder side of the material (bulk, far from the surface) to the hotter side (at the surface); upon arriving at the surface, the helium diffuses out of the tungsten at the free-surface boundary. The time needed to reach the steady periodic state is longer (∼3200 s for reaching 99% of the mean value of the oscillatory surface concentration) without accounting for the Soret effect than when accounting for it (∼1200 s). In both cases, due to the high thermal diffusivity of tungsten, both in its pristine and in its damaged state, the time needed for the species concentration to reach its steady periodic state is three orders of magnitude larger than the time needed for the temperature to achieve periodic steady state (∼2.1 s for reaching 99% of the mean value of the oscillatory surface temperature). The respective predictions for the helium content evolution when a constant equivalent heat load is used are 1% higher than those when the full ELM heat cycle is modeled. The discrepancy is attributed to the difference in the steady-state temperature profile: with the constant heat flux, the temperature depth profile is gradually changing, while this is not the case near the plasma-facing surface with the time-dependent heat load as can be seen in figure 2(c), and the Soret effect is sensitive to the local temperature gradient. Figure 4(b) shows the integrated deuterium content over the thickness of the tungsten slab as a function of time in a similar way as helium in figure 4(a). The same general remarks can be made about the deuterium behavior as for helium, with the difference that the deuterium content is higher than the helium content mainly because of the differences in their respective implantation profiles: for helium implantation, the peak of the implantation profile is located at a depth of 1.6 nm below the surface, while the respective peak is at 2.2 nm below the surface for deuterium [23]. The times to reach 99% of the periodic steady state are estimated to be ∼10 300 s without accounting for the Soret effect and ∼4400 s with the Soret effect accounted for, which are longer than those in the helium case because of the higher overall content for deuterium.

Tritium
Figure 4(c) shows the integrated tritium content over the thickness of the tungsten slab as a function of time, which is qualitatively similar with the helium content evolution, reported in figure 4(a). In a similar manner with the comparison of the helium with the deuterium content, the tritium content is higher than the helium content mainly because of the differences in their respective implantation profiles. The times to reach 99% of the steady periodic state are estimated to be ∼12 500 s without accounting for the Soret effect and ∼5400 s with the Soret effect accounted for, which are longer than those in the helium case because of the higher overall content for tritium. Despite both the hydrogenic species having the same implantation profile and heat of transport, the steadystate tritium content is a little higher than the deuterium content because of its smaller pre-exponential factor in its diffusion coefficient (see table 1).

Concentration profiles
The concentration profiles for helium and hydrogen (deuterium and tritium isotopes) as a function of depth in the tungsten slab are shown in figure 5. These profiles are obtained at the very end of the simulation, at the end of an ELM cycle. The results for the damaged tungsten case are not represented as we have established that the corresponding evolving temperature profiles do not vary much from those in the pristine tungsten case, resulting in practically indistinguishable species concentration profiles in the two cases. As observed previously, the concentrations of the hydrogenic species are higher than that of helium, with a more pronounced plateau over the depth of the material owing to the stronger temperature dependence of the diffusion coefficient (see table 1). The cases where the Soret effect is accounted for display a significantly lower concentration over the depth of the tungsten slab than in the absence of the Soret effect, although the concentration near the surface, where the impurities are implanted, are at the same level as when the Soret effect is not taken into account. The Soret effect seems to have a stronger impact for helium (due to its higher heat of transport than those of the hydrogenic species) with a sharper decrease of concentration (steeper concentration gradient) near the surface.
As mentioned before, the differences in the profiles between the equivalent heat load cases and full ELM heat model can be attributed to the fact that the temperature profile is linear at steady state with a constant heat flux. Therefore, at the steady-state condition under an equivalent heat load, the species flux balance is given by where γ is the fraction of implantation flux which is lost through the plasma-facing side and rear end of the tungsten slab, and β A is a positive constant, defined as Q * table 1). For simplification, we have assumed that species is implanted at a point coinciding with the peak depth and the entire domain is divided into two subdomains-from the plasma-facing side to the peak implantation depth and from the peak implantation deeper through the rest of the slab. We can conclude from equation (7) that the species concentration decays exponentially as a function of depth in the tungsten slab and the decay rate is proportional to the Soret coefficient, β A , consistent with the results in figure 5. For a complete analytical expression for the derived concentration profile see [26]. Finally, the predicted total content calculated by integrating the profile obtained using an equivalent heat load is ∼1% and ∼3% higher for helium and hydrogenic species, respectively, than when the full ELM heat cycle is modeled (see table 2) but leads to a drastic decrease in simulation time: for example, the time requirement to simulate 5000 s of tungsten response was reduced from 48 h to 20 min under identical computational conditions.
We have also computed the steady-state species fluxes diffusing out of the PFC material at the interface between tungsten and copper, at a depth of 6 mm from the plasma-exposed tungsten surface, without and with taking into account the Soret effect. These computed values are 2.6 × 10 14 m −2 s −1 and 6.1 × 10 12 m −2 s −1 , respectively, for helium, 7.1 × 10 13 m −2 s −1 and 3.7 × 10 12 m −2 s −1 , respectively, for deuterium, and 7.2 × 10 13 m −2 s −1 and 3.7 × 10 12 m −2 s −1 , respectively, for tritium, exhibiting a reduction factor of 43 and 19 for helium and hydrogen isotopes, respectively, when the Soret effect is included compared to the case where the Soret effect is not accounted for. The higher reduction factor observed for helium can be attributed to the larger value of its heat of transport compared to those of the hydrogen isotopes Table 2. A summary of periodic temperature and concentration profiles. min T and max T represent the minimum and maximum surface temperature (K) reached within an ELM cycle once the steady periodic state is established; similarly, min C and max C represent the total species content integrated over the depth of the tungsten slab (nm −2 ); τ ss,A is the time (s) required to reach 99% of the steady periodic state.  1). This reduction of the He and hydrogenic species fluxes at the interface between the PFC tungsten and the coolant pipe has major safety implications for fusion reactor design purposes. These results clearly show a pronounced impact of the Soret effect on the transport and retention of both helium and hydrogenic species, although it is important to note that the current simulations have not included the effect of selfclustering and trapping. Self-clustering of helium in tungsten is very strong and the literature contains many references showing that the presence of surfaces can induce trap mutation at lower helium cluster sizes than in the bulk [27][28][29][30]. The thermal-gradient-driven diffusion of helium towards the plasma-facing surface may influence the selfclustering and, therefore, the retention of helium. In addition, molecular-dynamics simulations and first-principles density functional theory calculations have shown a strong trapping interaction for hydrogen at the interface between a nanometersized, over-pressurized helium bubble and tungsten [31][32][33]. Furthermore, interactions between diffusing species and defects generated by 14 MeV neutrons also will affect species thermomigration. The use of atomistic modeling to evaluate the relative strengths of the interaction for self-clustering and trapping, versus the Soret effect on helium and hydrogenic species diffusional fluxes are required to further parameterize such behavior into our cluster-dynamics model. The results of such atomistic modeling will be reported in forthcoming publications, along with the extension of our clusterdynamics model and benchmarking the simulation results against experiments [34,35].

Summary
We have used NEMD simulation results from our previous work [5] to study the transport of helium and hydrogen in the presence of a thermal gradient in tungsten. Both helium and hydrogenic species have a negative heat of transport, which is implemented in our cluster-dynamics code Xolotl to analyze the consequences of the Soret effect on the retention of helium and hydrogen isotopes (deuterium and tritium) in pristine and damaged tungsten under ELM-like heat-load conditions. We found that, for all species examined, their total content over the depth of the tungsten slab simulated, after reaching their steady periodic state, was reduced (on average) by a factor of ∼3. While the actual value of the total content is strongly dependent on gas implantation and neutron irradiation conditions, especially on the location of the peak implantation profile [23], which would vary from one location to another in a divertor wall, the ratio of the species content between the cases with and without accounting for the Soret effect would be similar, under the assumption that transport of hydrogen isotopes is not heavily limited by self-clustering of helium and hydrogen binding to the resulting helium clusters. However, atomistic calculations have shown a strong propensity for helium self-clustering and overpressurized helium bubble formation, and a relatively strong trapping interaction of hydrogen with the helium clusters and bubbles. Therefore, the relative strengths of the interactions for self-clustering and trapping, versus the Soret effect on helium and hydrogenic species, would play a major role in determining the total species content under reactor-relevant conditions, which requires further investigation.
The use of an equivalent heat load showed encouraging results for justified simplifications in future simulations with drastic improvements in terms of computational cost reduction. Finally, the results with near-surface damage in tungsten did not display any distinguishable differences with respect to those when pristine tungsten was considered, since we only included the effect of near-surface damage on thermal properties. In the future, we plan to extend the simulations to take into account species clustering and trapping, with the possible addition of pre-existing damage to model the nearsurface damaged tungsten more accurately. For those future simulation efforts, use of an equivalent heat load in place of an ELM-like periodic heat load would dramatically reduce the corresponding computational costs.