1. Introduction
Accurately estimating the quantity of tritium retained in plasma facing components (PFCs) is a key safety requirement for licensing future fusion power reactors, well as for the International Thermonuclear as Experimental Reactor (ITER) [1].  These components not only see the very harsh environment produced by a plasma, such as surface erosion, severe surface heat flux, and high rates of plasma ionic specie implantation, but the 14 MeV neutrons generated by the fusion process create significant damage to the PFC material lattice by atomic displacements and elemental changes through transmutations [2].  This damage acts as ideal traps for tritium retention as implanted tritium diffuses towards the surfaces of the PFC, where it is either released back into the vacuum vessel, or into the PFC’s coolant. Tritium retention in PFC materials is an area of experimental research at the Idaho National Laboratory (INL) under the US/Japan TITAN collaboration.  Recent experiments with the Tritium Plasma Experiment (TPE), located in the Safety and Tritium Applied Research (STAR) facility at the INL, demonstrate that this damage can only be simulated by codes like the Tritium Migration Analysis Program (TMAP) [4] if multiple trap sites are used to simulate the damage generated by The TMAP code has been released as neutrons [3]. various versions, the most recent of which has the capability of modeling at most three user specified trap site types. However, even with this number of sites, TMAP was unable to adequately simulate deuterium release from tungsten targets exposed to TPE plasmas that had first been irradiated to a damage level of 0.025 dpa in the High Flux Isotope Reactor (HFIR) at the Oak Ridge National Laboratory (ORNL). In this paper we describe recent extensions of the TMAP trap site model to treat as many traps as required by the code to simulate the retention of tritium in neutron damaged tungsten material.  In the sections that follow, we describe the equations for these new trap site models, the experiments conducted at t the INL to develop the data needed to benchmark the implementation of these models, and the application of these models to deuterium retained in neutron irradiated tungsten after exposure to a TPE plasma.
2. Trapping Models in TMAP
The TMAP calculates the time-dependent response any number of gaseous species in a system of solid structures or walls, and related gas filled enclosures by including: . . Implantation of species in structures, of species through structure Movement surfaces, governed by molecular dissociation and atomic recombination, or by a solution law such as Sieverts’ or Henry’s Laws . . . Movement in the structure as treated by Fick's law of bulk diffusion with the possibility of specie trapping in material defects, Thermal response of structures to applied heating or boundary temperatures, Chemical reactions within the enclosures, and . specified convective between User flow enclosures Because the equations used to simulate these processes are highly nonlinear, a Newton-Krylov type solver is used to update these equations at each point in time during a calculation. The basic equation for movement of these species in the bulk material is the one-dimensional Fick's Second Law of diffusion written as follows:
where: . C m is the mobile specie concentration (m 3 ) .
. n is the number of traps allowed by TMAP The rate at which the mobile species are entering trap defects in the bulk material is: where: α , is the trapping rate coefficient (s 1 ) . . .
The definition of these terms i (3) and λ is the jump distance or lattice constant (m), c ,  is the trap site concentration, N is the bulk material number density (atoms-m 3 ), v o is the Debye frequency (s 1 ), E t is the trap binding energy (eV) and k is Boltzmann’s constant (eV/K). In order to allow the use of more than the three traps presently allowed by TMAP7, it was a trivial matter of merely changing the TMAP code to accept more than traps, which was also done for the purpose of this study. However, when viewing the Thermal Desorption System (TDS) data from neutron irradiated tungsten it became apparent that perhaps a more accurate approach would be to modify the TMAP code trapping model to accept any distribution of trap site energies; in particular for this study, a distribution based on previous work performed by the INL Fusion Safety Program (FSP) in developing aerosol resuspension models, which is a lognormal or bimodal lognormal distribution [5]. The density function for a lognormal trap distribution in terms of trap energy is written as:
where μ is the lognormal mean and σ is the lognormal standard  deviation. Given this distribution, the total concentration of trap sites (c ? ) at any location in the bulk material equals:
and the trapped concentration of a given specie C t at that location, and at a given point in time equals:
Differentiating 6  with respect Eq. to time and substituting from Eqs. 2 and 3, the result for the rate of change of the total trapped concentration of a given specie becomes:
While Eq. 7 could be used directly in Eq. 1, the TMAP code has always been meant to be a relatively fast running computer code capable of performing quick safety scoping and data reduction analysis tasks for the INL FSP.  To accomplish this for the present trapping model, we again adopted the same approximation made for our dust resuspension model.  By noting that high energy traps will be the last traps to release their trapped inventory, we assume that the traps of this distribution will fill from high energy to low energy in that order, that is from i = N to 1 as illustrated in Fig. 1.  The advantage of this assumption is the TMAP is now only tracking the total trapped inventory in this distribution. Since the fraction of traps filled at any point in time equals the trapped inventory divided by the trap density (e.g., F = C t/c t ), finding the trap site that is changing in concentration for this approximation is a simple matter Cumulative of interpolating  with respect to the Distribution Function (CDF) for lognormal a distribution. This trap is assigned the net transient trapping and release rates for the entire distribution as its concentration is changing in time.
3. Experiments Conducted at t the INL
Under the US/Japan TITAN collaborative agreement, tungsten samples that are 6 mm in diameter and 0.2 mm thick were irradiated to three different damage levels (0.025, 0.3 and 2.4 dpa) in the HFIR reactor at ORNL at a temperature of 50°C. Following an appropriate cool down period, these samples were exposed to the plasma of the Tritium Plasma Experiment (TPE) at INL's STAR facility.  The TPE is a linear plasma device [6] capable of exposing target material to low-energy (100 eV) deuterons at a plasma ion flux of up to (5–7)×10 1  m-² s-¹.  For the (2) irradiated tungsten samples the plasma was maintained until an ion fluence of (5–7)×10 25 m-2 was achieved for sample temperatures of 100, 200 and 500°C. After the deuterium plasma exposure, the samples were removed from TPE and were transferred to the TPE Thermal Desorption System (TDS) in air, where the samples were heated to produce a controlled linear temperature change from 25 to 800°C.  The off-gas from the samples was passed through the Residual Gas Analyzer (RGA) of the TDS that was specially calibrated to measure the absolute concentration of off-gassing D 2 . Only the peaks associated with HD (mass 3) and D 2 (mass 4) were considered in this analysis. Fig. 2 presents a TDS spectrum obtained for a 0.025 dpa sample used in this study. Several of the exposed samples were sent to the Ion Beam Laboratory at University of Wisconsin, Madison for Nuclear Reaction Analysis (NRA) [7] prior to undergoing assay in the STAR facility’s TDS.  At this facility, a high energy 20 mA helium ion beam is directed at t the target to create a D( 3 He, p) 4 He nuclear reaction with the deuterium implanted by a TPE plasma. The high energy protons emanating from the sample are used to estimate the depth profile of the deuterium retained in the sample.  Given the relatively immobile nature of deuterium in tungsten at room temperature, this depth profile can also be used to infer the hydrogen trap site concentrations for those traps that saturate under TPE irradiation conditions.  Fig. 3 contains a sample of NRA data as used in this study.
4. Application of Trapping Models to Irradiated
Tungsten Deuterium Retention Data
component examination A of any of key experimental data is the validation of the physics models being solved by the computer code being applied. In As was the case for Ref. [3], in this study we adopted the hydrogen diffusivity (corrected for deuterium) and solubility formulas proposed by Frauenfelder [8], where (D = 2.9×10 7 exp(-0.39 eV/kT) (m /s) and K = 1.79×10 24 exp(-1.04 eV/kT) (D/m 3 -Pa 1/2 ).  In literature, surface recombination coefficient formulas can vary by several orders of magnitude [9].  However at this time, the only formula based on experimental data for tungsten is that by Anderl of K, = 3.2×10 15 exp(-1.16 eV/kT) order to verify TMAP’s trapping models against TPE results, property data for tungsten diffusivity, solubility, and surface recombination are required. (m /s) was adopted for this study as a representative value for a plasma vacuum environment [10]. Physical data of trap density, D + implantation depth, plasma surface flux intensity and target temperatures are also required.  The NRA measured D/W atom fraction presented in Fig. 3 was used as a starting value for the total trap site density in neutron irradiated tungsten to a damage level of 0.025 dpa.  An implantation depth of 3 nm, as estimated by Wampler [2], was used in our study.  In addition, in order to understand if these new TMAP models accurately address the physics of plasma driven retention in tungsten, the actual TPE plasma operating history for the TDS data being examined must For the data presented in Fig. 2, the TPE be used. plasma flux (ramped up to 5×10 21 D*/m²-s and held for 7200 s) and target temperature (under plasma heating to 200°C) histories for the TPE run that resulted in the TDS data of Fig. 2 were follow in the TMAP model as reported by Shimada [3]. Our first goal is to verify the TMAP results obtained by Shimada [3], for this same TDS case. A good comparison between TMAP7 and TDS data was obtained by these authors by assuming six traps. However, because that version of TMAP only allowed a maximum of three user specified trap sites, the results of two separate TMAP simulations had to be manually added to obtain the result presented.  The trap energies and total density for these traps were reported as 0.9, 1.1, 1.3, 1.5, 1.75 and 2.0 eV, and 1.3 at% D/W, respectively. As noted by these authors, the trap density required to match the TDS was much higher than the NRA D/W atom fraction measured for the 0.025 dpa sample that underwent the NRA in Fig. 3.  There are a number of possible explanations for this difference, including losses of deuterium by diffusion and surface release between the time of TPE irradiation and NRA examination, uncertainty in the NRA results, or differences in TPE operating history for the TDS sample of Fig. 2 and NRA sample of Fig. 3, etc.  However, because we do not have the NRA depth profile for the sample used in the TDS measurement of Fig. 3, it is equally plausible that the trap density is 1.3 at% for the sample analyzed by the TDS. Fig. 4 contains the results of using six traps used in this study to simulate the TDS of Fig. 3 using the new multi-trap site TMAP model.  Unlike the simulations of Ref. [3], where the decision to arbitrarily assign the trapped deuterium concentration profile as the initial condition for performing their TDS simulations, this result was obtained by simulating the implantation process of the TPE plasma, with the goal of determining if the adopted trap sites would actually saturate, as assumed by Ref. 3, at the irradiation temperature (200°C) The six traps required to experienced by the target. simulate the TDS spectrum had energies of 0.9, 1.0, 1.1, 1.2, 1.4 and 1.6 eV, and a total trap density of 0.35 at%. One aspect of this simulation that was not true to TPE plasma operation details deals with the filling of the lower energy trap sites.  TMAP predicts that at 200°C the 0.9 eV traps will not saturate.  Since TMAP predicts that the mobile deuterium concentration in the tungsten target, which is several orders in magnitude less than the trapped concentration, is too low to fill these traps once the TPE plasma has been terminated, the TPE plasma ramp down had to be extended in time beyond what had actually occurred.  The result marked ''no overfill'' in Fig. 4 illustrates the predicted TDS spectrum when the plasma was ramped down as reported.  While additional validation of this new TMAP capability is still needed to verify the accuracy of this result, it does suggest that the deuterium being released in the TDS spectrum below 200°C had nothing to do with trapping on dislocations but is more likely associated with the deuterium that resides in surface blisters known to develop in tungsten at the implantation fluence experienced by this TPE target (note discussion for Fig. 4 of Ref. [2]). Fig. 5 presents the results obtained for this same case with the new TMAP lognormal trap model.  An adequate representation of the TPE TDS spectrum was obtained with a truncated (at 1.6 eV) bimodal lognormal distribution, which is a distribution that is the sum of two lognormal populations with proportions p and q, namely:
φ(E)= pφ)(E:μι,σι)+ qφ2(E:μ2,σ2)
(8) For the comparison in Fig. 5, p and q were both set to 0.5, and the mean and standard deviations of the lognormal distribution φ 1 were set at 0.95 eV and 0.18 eV and for φ 2 set at 1.45 eV and 0.21 eV.  The total trap density for this case was 0.47 at%.  As can be seen, the simulation represents the data well.  As was the case for the six trap model, the lower energy traps could only be filled by ramping TPE down more slowly that what actually occurred.
5. Summary and Conclusions
Accurately estimating the quantity of tritium retained in or permeating through PFCs is a key safety requirement for licensing future fusion power reactors. The magnitude of these processes is greatly affected by damage created in the structural materials of these PFCs by the 14 MeV neutrons generated by the fusion process. This damage acts as ideal trap sites for tritium retention within the PFCs. In this paper, we have introduced two extensions to the TMAP computer code trap site model. These extensions allow the TMAP user to specify as many trap sites or a distribution of trap sites as needed to simulate the material damage being modeled.  These models were benchmarked against data generated under the US/Japan TITAN collaborative research agreement in the INL's TPE device for retention of deuterium in tungsten that had been irradiated to a damage level of 0.025 dpa in the ORNL HFIR. We found that the multi-trap extension to TMAP could reproduce the off-gassing data from a tungsten sample exposed to a TPE plasma by employing six traps of energies of 0.9, 1.0, 1.1, 1.2, 1.4 and 1.6 eV, and a total trap density of 0.35 at%.  Alternatively, the data could also be matched if a bimodal lognormal trap site distribution is used that has proportions p and q equal 0.5, and mean and standard deviations of the lognormal distribution φ 1 of 0.95 eV and 0.18 eV and for φ 2 of 1.45 eV and 0.21 eV, respectively.  The total trap density for this case was 0.47 at%. While additional validation of these new TMAP capability is still needed to verify the accuracy of this result, the results appear to suggest that the deuterium being released in the TDS spectrum below a temperature of 200°C had nothing to do with trapping on dislocations but is more likely associated with the deuterium that resides in surface blisters known to develop in tungsten at implantation fluence experienced by this TPE target. Future work in this area with a broader set of experimental data must be undertaken to confirm this conclusion.
Acknowledgments
This work was prepared under the auspices of the US Department of Energy, Office of Fusion Energy Science, under the DOE Idaho Field Office contract number DE- AC07-05ID14517.