VALIDATION OF THERMAL SCATTERING LAWS FOR LIGHT WATER AT ELEVATED TEMPERATURES WITH DIFFUSION EXPERIMENTS

Light water is the most important neutron moderator in many reactor applications. Thermal neutron scattering kernels for H bound in H 2 O are represented by ENDF-format thermal scattering laws (TSLs) evaluated at discrete temperatures T . Nuclear data evaluations are conventionally validated utilizing a combination of critical benchmarks and experimental cross section data. Existing public critical benchmarks may be inadequate for testing water TSLs over the full range of T of interest in reactor applications, and experimental thermal scattering cross section data for water is sparse at elevated T . In this work, MC21 is used to simulate the decay of the equilibrium thermal neutron flux in light-water spheres of several radii. The fundamental-mode time eigenvalue α is calculated at 22 °C and 227 °C (at saturation pressure) for each sphere using three different H-H 2 O TSLs. Fitting α to a polynomial function of geometric buckling allows calculation of the thermal neutron diffusion length L . Extrapolation lengths are treated as a function of the transport mean free path and geometry. Experimental L data from 25 publications at 49 temperatures (from 10 °C to 295 °C), computed by several different time-dependent and space-dependent decay methods, is used to develop an empirical fit of L vs. T . The MC21-calculated L for the TSLs tested are within ±1% of the predicted values at 22 °C and 227 °C. This validation approach, which may be repeated at arbitrary T , constitutes an integral benchmark specific to and characterizing the detailed physics of the thermal scattering kernel applied. of nuclear power. Accurate modeling of neutron flux distributions and neutron multiplication in water-moderated thermal systems requires thermal scattering cross sections that account for the effect of chemical binding and interactions of H 2 O molecules. Early efforts to characterize this effect included extensive measurements of the diffusion properties of neutrons in water. Consequently, there is a large quantity of experimental data available regarding thermal neutron diffusion in water from several dozen publications over the large range of temperatures of interest in reactor applications.


INTRODUCTION
Light water is the most common neutron moderator utilized in a wide range of reactor designs and criticality safety applications, and thermal neutron scattering in water has been heavily studied since the earliest days of nuclear power. Accurate modeling of neutron flux distributions and neutron multiplication in watermoderated thermal systems requires thermal scattering cross sections that account for the effect of chemical binding and interactions of H2O molecules. Early efforts to characterize this effect included extensive measurements of the diffusion properties of neutrons in water. Consequently, there is a large quantity of experimental data available regarding thermal neutron diffusion in water from several dozen publications over the large range of temperatures of interest in reactor applications.
Integral and differential thermal scattering cross sections can be computed from the thermal scattering law (TSL) characterizing the interatomic structure and dynamics of a material in terms of neutron momentum and energy transfer probability. The TSL is conventionally represented in ENDF-format [1] libraries by the tabulated ܵ(ߙ, ߚ, ܶ) function, where ߙ is a momentum transfer factor, ߚ is an energy transfer factor, and ܶ is material temperature. Several TSLs for H bound in light water (H-H2O) have been developed using progressively improved methods. The early 1968 General Atomics (GA) model by Koppel and Houston [2] was based on differential scattering measurements by Haywood and Thorson at 20 °C and 150 °C in 1962 [3]. The 1984 IKE model by Keinert and Mattes [4] was based on differential scattering measurements by Page and Haywood at 295 K and 550 K in 1968 [5]. It included a more detailed ߙ and ߚ mesh and showed better consistency with a variety of experimental cross section data than the GA model. In the United States ENDF nuclear data library, various versions of the GA model were used in ENDF/B-III (released in 1972) through ENDF/B-VI. 8 (released in 2001). The GA model was then replaced with an improved version of the IKE model [6] in ENDF/B-VII.0 [7] (released in 2006). In 2017, ENDF/B-VIII.0 [8] was released using the CAB model by Márquez Damián et al. [9]. The CAB model uses molecular dynamics (MD) simulations to generate the TSL and shows further improvements against experimental cross section data, particularly for cold neutrons (< 1 meV) in water at room temperature (RT).
The ܵ(ߙ, ߚ, ܶ) function is tabulated on a discrete grid of 18 temperatures in the CAB model (ranging from below RT to above the supercritical temperature for water) and on a grid of 9 temperatures or fewer in earlier models. For every ܶ, an underlying TSL physics model must be separately defined that depends on interpolation/extrapolation of experimental measurements and/or on tuning/adjustment of MD force field parameters. Nuclear data evaluations are conventionally validated using a suite of ICSBEP integral critical benchmarks [10] along with experimental cross section data. IRPhEP benchmarks [11] have also been recently used. The vast majority of ICSBEP and IRPhEP benchmarks are at RT only, and the exceptions (LEU-COMP-THERM-026/032/046, KRITZ-LWR-RESR-001/002/003, etc.) include limited data at elevated ܶ. Likewise, experimental thermal scattering data for water at elevated ܶ is very sparse. The little data available is either for cold neutrons (< 1 meV) [12,13,14] of marginal interest in reactor applications or for limited select ܶ [15,16,17]. Not surprisingly, light-water TSL performance has been thoroughly validated at RT, but comprehensive non-proprietary validation across the full range of ܶ of interest in reactor applications is lacking. Moreover, due to the interpolated/extrapolated ܶ-dependent physics models, there is no assurance that good performance at one ܶ implies good performance at another.
Notably, a significant non-physical ܶ-dependent reactivity bias was identified by Naval Nuclear Laboratory when modeling the Neptune experiment [18] from 27 °C to 60 °C with an early beta version of the CAB light-water TSL. This bias was not observed with other TSL models. While Neptune model details are proprietary, the Neptune results were presented at the 2017 OECD-NEA Working Party on International Nuclear Data Evaluation Cooperation (WPEC) meeting [19]. CAB subsequently modified the TSL physics models, and retesting indicated the Neptune bias was mostly eliminated.
The neutron diffusion length (or relaxation length) ‫ܮ‬ is an integral property of the neutron scattering and absorption characteristics of a material. Comparing ‫ܮ‬ calculated with H-H2O TSLs to the substantial experimental data for ‫ܮ‬ in light water available over a wide range of ܶ, independent of the influence of other materials, may be better suited to validating the TSL physics models than critical benchmarks.

THEORETICAL AND EXPERIMENTAL ASPECTS OF NEUTRON DIFFUSION
Many types of experiments have been conducted to measure thermal neutron diffusion properties in water. The two primary classes of diffusion experiments are stationary (space-dependent) methods and nonstationary (time-dependent) methods [20,21]. For a neutron flux in thermal equilibrium, both methods can be used to determine ‫,ܮ‬ which is related to the average straight-line distance from a source (or starting point) to the point of absorption as neutrons diffuse in a homogeneous medium.
For the non-stationary case, the time-dependent one-speed neutron diffusion equation for a homogeneous medium is The variables ‫,ݒ‬ ߮, , ‫,ݐ‬ ‫,ܦ‬ ߑ ୟ , and ܵ are defined as follows: neutron speed, neutron flux, spatial position, time, diffusion coefficient, macroscopic absorption cross section, source strength. For a non-multiplying system with no neutron source, the fundamental-mode solution for Eq. (1) has the form where the time eigenvalue ߙ is given by The term ‫ܤ‬ ଶ is defined as geometric buckling and is simply the fundamental-mode eigenvalue ‫ܤ‬ ଵ ଶ .
Geometric buckling characterizes the effect of boundary leakage on the flux spatial curvature and on the flux decay rate in time. For spherical systems of radius ‫ݎ‬ , where ‫ݖ‬ is the extrapolation length, or the distance at which a tangential extrapolation of the diffusion-theory spatial flux profile at a vacuum boundary goes to zero.
It is known that diffusion theory does not perform well in small systems, near boundaries, near strong sources or absorbers, or for rapid transients on the order of the mean collision time. Moreover, the onespeed model is a very coarse approximation. Nevertheless, for a real homogeneous source-free nonmultiplying system with a thermalized neutron flux, Eq. (2) still describes the time evolution of the fundamental-mode flux profile. For a physical case governed by Eq. (2), the time eigenvalue ߙ can be expanded as a Taylor series in ‫ܤ‬ ଶ [22,23,24] in terms of diffusion properties: In this formulation, ܽ ୣ is an "effective" absorption term, ‫ܦ‬ ୣ is an "effective" diffusion coefficient, and the one-speed approximation is relaxed (ܽ ୣ and ‫ܦ‬ ୣ are spectrally averaged). The parameter ‫ܥ‬ is a diffusion cooling coefficient associated with preferential "transport" and boundary leakage of higher-energy neutrons in the thermal spectrum. The coefficients in buckling are functions of the integral and differential thermal neutron scattering properties of the material. The second-order and higher terms in buckling are essentially "transport corrections" to the diffusion-theory time-eigenvalue form given by Eq. (3).
A common method for determining these physical diffusion parameters is the pulsed-neutron die-away (PNDA) experiment [20,21]. A material sample is bombarded with short (on the order of microseconds) high-energy neutron pulses from a neutron generator. After each pulse, the system neutron flux magnitude is measured at a detector position as a function of time. After sufficient waiting time, the neutron flux is thermalized, higher harmonics are negligible, and the asymptotic ߙ can be evaluated. By determining ߙ over a range of ‫ܤ‬ ଶ (or for samples of varying dimension), the parameters ܽ ୣ , ‫ܦ‬ ୣ , and ‫ܥ‬ can be calculated with a quadratic fit of ߙ to ‫ܤ‬ ଶ (or higher-order fits may be applied as appropriate). PNDA experiments have been used extensively to study the diffusion properties of water [23,[25][26][27][28][29][30][31][32][33][34]. If an external source makes up for losses in the flux such that it is constant in time, this is a stationary case. The relevant one-speed diffusion equation for a homogeneous medium is where ‫ܮ‬ ≡ ඥ‫ߑ/ܦ‬ ୟ . For an infinite homogeneous non-multiplying system with a point source at ‫ݎ‬ = 0, the fundamental-mode flux profile solution to Eq. (5) converges as ‫ݎ‬ → ∞ to the form where ߢ = ‫.ܮ/1‬ For a real system, Eq. (6) still describes the asymptotic solution far from a source or boundary, and the spatial eigenvalue ߢ will depend on the physical ‫.ܮ‬ Experimentally, one common method of determining ‫ܮ‬ for this stationary case is by using a 25 keV Sb-Be photoneutron source and indium activation foils to measure the spatial flux decay, with appropriate source corrections, far from a boundary [35,36]. Other stationary experiment techniques exist as well, including the poisoning technique commonly used to evaluate the spectrally averaged ߑ ୟ for the pure medium [36,37].
In Eq. (3), it can be seen that a solution of ߙ = 0 requires ‫ܤ‬ ଶ = ‫ܮ/1−‬ ଶ , since ‫,ݒ‬ ߑ ୟ , and ‫ܦ‬ are physical properties that cannot be negative. Recalling that ‫ܤ‬ ଶ characterizes boundary leakage, this situation is equivalent to stating that net leakage is negative due to an external source. However, we have noted that Eq. (4) provides a more physical representation of ߙ. In this case, ‫ܤ‬ ଶ = ‫ܮ/1−‬ ଶ still gives the required negative net leakage to yield ߙ = 0 with a fixed source present [23,36]. Solving the quadratic expansion of Eq. (4) in the stationary case by substituting ‫ܮ/1−‬ ଶ for ‫ܤ‬ ଶ gives While this solution for ‫ܮ‬ is based on the stationary case, ‫ܮ‬ is a material property. Therefore, Eq. (7) gives a more physical form for ‫,ܮ‬ in terms of expanded diffusion coefficients, that applies in non-stationary and stationary cases. When 4ܽ ୣ ‫ܦ/ܥ‬ ୣ ଶ ≪ 1 (which is true for water), Eq. (7) may be closely approximated by the more familiar form ‫ܮ‬ = ඥ‫ܦ‬ ୣ /ܽ ୣ [23,34]. Indeed, Honeck [38] has shown that the direct determination of ‫ܮ‬ from spatial flux decay measurements is theoretically equivalent to the determination of ‫ܮ‬ from expanded diffusion coefficients calculated by temporal flux decay measurements. Koppel et al. [39] and Parks et al. [40] demonstrate this equivalence by comparing experimental results from both classes.
Of the two classes of diffusion experiments, the stationary method allows more direct access to ‫,ܮ‬ although geometric flux corrections are required. If desired, the expanded diffusion coefficients given in Eq. (4) can be equivalently determined in the stationary case using the poisoning technique [41]. The non-stationary PNDA method allows more direct access to the expanded diffusion coefficients but can be subject to ambiguities in fitting and in the calculation of ‫ݖ‬ and ‫ܤ‬ ଶ [35]. Both methods have been extensively used.

Temperature Dependence of the Diffusion Length in Light Water
The temperature dependence of ‫ܮ‬ in water has been studied extensively by various stationary and nonstationary experimental methods [28,30,31,[35][36][37][42][43][44], and several different empirical functions have been proposed [28,31,35,43,44]. The physical variation of ‫ܮ‬ with ܶ depends on the non-linear relationship of density (ߩ) vs. ܶ as well as on changes in the dynamics of H2O molecules (such as changes in clustering and rotational hindrance). Under the assumption that microscopic cross sections are ܶ-independent and have ‫ݒ/1‬ dependence over a Maxwell-Boltzmann neutron spectrum, diffusion theory predicts that ‫ܮߩ‬ should vary as √ܶ [42,44]. Physically, experimental data has been shown to follow the form ‫ܮߩ‬ = ܶ ோ , where ܴ is slightly less than 0.5. Reier and de Juren [35] give ܴ = 0.4343 for measurements from 23 °C to 244 °C, while Besant and Grant [44] give ܴ = 0.487 for measurements over the narrower range of 24 °C to 82 °C. Each of these empirical values was derived from only one experimental data set.
Not surprisingly, clear systematic biases exist between different publications, where the data spread exceeds what may be reasonably expected based on the authors' quoted uncertainties. Moreover, ‫ܮ‬ data from the 10 PNDA sources in Fig. 1 falls (on average) 0.8% below the empirical curve, while data from the 15 stationary method sources falls (on average) 0.2% above (and has half the average quoted uncertainty).
Taking into account the ambiguity of the observed systematic biases, a best estimate of 0.5% is given for the uncertainty in calculating ‫ܮ‬ with the empirical function ‫ܮߩ‬ = (0.1766 g/cm ଶ )ܶ .ସ଼ଷ (for ܶ in Kelvin).

Physical Treatment of the Extrapolation Length
The procedure for computing diffusion coefficients and ‫ܮ‬ for the non-stationary case has been laid out in Section 2. However, results will evidently depend on treatment of the extrapolation length ‫ݖ‬ , which is an artifact of diffusion theory. For a system of arbitrary shape and dimensions, there is always an "effective"  ‫ܤ‬ ଶ that will satisfy the solution to the fundamental-mode time eigenvalue ߙ given by Eq. (4). Therefore, the "correct" value for ‫ܤ‬ ଶ is the value that reproduces the physical ߙ in a real homogeneous system with material properties ܽ ୣ , ‫ܦ‬ ୣ , ‫,ܥ‬ etc. Likewise, the "correct" ‫ݖ‬ in any standard-geometry ‫ܤ‬ ଶ formula is that which yields the "correct" ‫ܤ‬ ଶ [54].
The treatment of ‫ݖ‬ and ‫ܤ‬ ଶ in PNDA experiments have been the subject of considerable study [24,[54][55][56][57][58], including for light water and for spherical systems. The appropriate ‫ݖ‬ is a function of material density, material scattering and absorption properties, and geometry. This relationship can be expressed as ‫ݖ‬ = ߣ ୲୰ ‫ܨ‬ . The term ߣ ୲୰ is the transport mean free path. The geometry factor ‫ܨ‬ is shape dependent and is a nearly-constant function of ‫ܤ‬ ଶ , except for small geometries less than a few mean free path lengths (ߣ ୲୰ ≈ 0.5 cm in water at RT). In the present work, the ‫ݖ‬ values computed by Schmidt and Gelbard [55] for lightwater spheres at RT (~22 °C) are used to calculate ‫ܤ‬ ଶ , as they appear to be the most accurate for the range of buckling examined (including for smaller radii where corrections for flux curvature become significant).
To determine appropriate ‫ݖ‬ at 500 K (the elevated temperature tested in this work), we note that the theoretical value for ‫ܮ‬ ହ is given by (0.1766 g/cm ଶ )(500

RESULTS AND DISCUSSION
To validate the performance of ENDF light-water TSLs, PNDA simulations were conducted with the MC21 Monte Carlo transport code [59] at 22 °C and 227 °C (at saturation pressure) for eight spheres of water with radii ranging from 3.6 cm to 14.4 cm. (295.15 K and 500 K exactly were used.) This method was selected because PNDA simulations with MC21 were previously found to be highly successful in validating a theoretically-generated TSL for ice Ih over a wide temperature range [60]. Three ENDF light-water TSL models for H bound in H2O (H-H2O) are tested in this work: ENDF/B-VII.1 [61], ENDF/B-VIII.0 [8], and a version under development by North Carolina State University (NCSU) based on MD simulations [62,63].
Oxygen is treated as a free gas in all cases. The spherical radii, extrapolations lengths, geometric buckling, and calculated ߙ eigenvalues at 22 °C and 227 °C for the three ENDF TSLs are given in Table I. Eigenvalue results treating hydrogen as a free gas are also given at 22 °C. The calculated results for ‫ܮ‬ and the expanded diffusion coefficients (up to second order) are given in Table II and compared to experimental averages. For the free gas case at 22 °C, the calculated ‫ܮ‬ is 3.603 cm (30% higher than the experimental average). A graphical display of the calculated ‫ܮ‬ results (and experimental data) in the context of the empirical ‫ܮߩ‬ = (0.1766 g/cm ଶ )ܶ .ସ଼ଷ fit is shown in Fig. 2. The experimental data points at exactly 22 °C were removed in Fig. 2 to allow easier visualization of the MC21 results.

Uncertainties
Many factors can influence the uncertainty in the calculated values of diffusion parameters with Eqs. (4) and (7). These include the order and method of the fit, the number of ߙ values and range of ‫ܤ‬ ଶ used, the values of ‫ݖ‬ used to calculate ‫ܤ‬ ଶ , statistical uncertainty in ߙ, room return (for experiments), irregularities in system shape (for experiments), the treatment of oxygen as a free gas (for simulations), and others. Quantifying the effect of these factors can be difficult. It is also evident from examining many independent sets of experimental data (and quoted uncertainties) that significant unidentified systematic biases are common. Of course, this applies to stationary methods as well.   The experimental values and uncertainties for ܽ ୣ , ‫ܦ‬ ୣ , and ‫ܥ‬ are unweighted and based on various combinations of data at or near RT from Refs. 23, 25-27, 29-32, 39, and 41. All data for ܽ ୣ and ‫ܦ‬ ୣ is adjusted to 22 °C. Data for ‫ܥ‬ is unadjusted. The experimental value and uncertainty for ‫ܮ‬ at 227 °C is based on the empirical fit from Section 2.1. The estimated uncertainty in MC21-calculated using the given PNDA method is 0.3%. In Table II, with the exception of ‫ܮ‬ ୰ୣ , diffusion parameter experimental averages are calculated as a straight average from numerous publications without considering the quoted experimental uncertainties (due to obvious systematic biases). Likewise, uncertainties given for these averages are based solely on the statistical spread of available data. Qualitatively, the experimental average for ‫ܮ‬ ୰ୣ at 22 °C is likely very accurate for two reasons. First, a large number and variety of sources are used in its computation. Second, ‫ܮ‬ is an integral quantity and is quite insensitive to higher-order diffusion effects for light water.
Uncertainties in the average experimental values for ܽ ୣ and ‫ܦ‬ ୣ are considerably higher than for ‫.ܮ‬ This is partly because ‫ܮ‬ ≈ ඥ‫ܦ‬ ୣ /ܽ ୣ , but also because fewer publications had this data available. Furthermore, much of the experimental ‫ܮ‬ data was directly calculated with static methods and did not result from a polynomial fit. The experimental ܽ ୣ , ‫ܦ‬ ୣ , and ‫ܥ‬ data in Table II is primarily from PNDA sources. Considering that the PNDA-based experimental ‫ܮ‬ data used in the work is biased about 1% lower than the stationary-method ‫ܮ‬ data (see Section 2.1), it is possible that the experimental ‫ܦ‬ ୣ average is biased too low (since ‫ܮ‬ ≈ ඥ‫ܦ‬ ୣ /ܽ ୣ and ܽ ୣ can be determined more accurately than ‫ܦ‬ ୣ ). The uncertainty in experimentally calculated ‫ܥ‬ (and in higher-order terms) is always large [23,29,30,33,34]. This is because higher-order diffusion effects are very weak in light water (i.e., the absorption cross section is very small compared to the scattering cross section, and differential scattering probabilities have near-linear dependence on the scattering angle). Therefore, effectively probing ‫ܥ‬ requires the use of extremely small geometries (or extremely high buckling), which can be impractical.
Due to the simplicity of the system models, it is relatively inexpensive to calculate ߙ in MC21 with very small statistical uncertainty. For the values given in Table I accuracy of the ߙ fit to buckling, a second-order fit was first performed over several combinations of lowbuckling cases, where the influence of higher-order terms in ‫ܤ‬ ଶ should be nearly zero. This allowed a very accurate determination of ܽ ୣ , or the ‫-ݕ‬intercept. Next, to determine ‫ܦ‬ ୣ and ‫,ܥ‬ a second-order fit was performed over all data with ܽ ୣ held constant. Higher-order fits were considered but provided no benefit. Higher-order fits would be necessary if a larger range of ‫ܤ‬ ଶ is used.
It should be noted that the ENDF/B-VII.1 covariance data gives an uncertainty of 2.5% on thermal ߪ ୟ ( 1 H) at 0.0253 eV [61], while ENDF/B-VIII.0 gives 2.1% [8]. These uncertainties are likely overestimated. In ENDF/B-VIII.0, ߪ ୟ ( 1 H) = 0.3327 barns. This is essentially identical to the value given in Ref. 64 Based on the known statistical uncertainties in ߙ, assessing the sensitivity to and influence of likely errors in ‫ݖ‬ , examining various fitting strategies, and assuming the uncertainty in thermal ߪ ୟ ( 1 H) is 0.2%, it is estimated that the uncertainty in the MC21-calculated ‫ܮ‬ in this work (due to these factors) is 0.3%. The uncertainties in calculating the expanded diffusion coefficients are not considered but will be larger. Any statistical uncertainties in calculating diffusion parameters could be arbitrarily reduced by simply increasing neutron histories and calculating more ߙ points over a larger range of ‫ܤ‬ ଶ . No consideration is made regarding the impact of treating oxygen as a free gas. However, given the results in this work (and considering its very small cross section and high mass relative to 1 H), the bias is likely extremely small.

CONCLUSIONS
The neutron diffusion length ‫ܮ‬ is an integral physical property of the scattering and absorption characteristics of a material. In this work, the ability to accurately calculate ‫ܮ‬ for thermal neutrons in light water with MC21 PNDA simulations at very different temperatures has been established by leveraging a large quantity of experimental data from numerous publications using a variety of experimental methods. materials, the models are extremely simple and fast to run, and a large quantity of experimental data for thermal neutron diffusion in light water is available for comparison. Repeating the given method over a fine mesh of temperatures should result in a smooth inflectionless curve of ‫ܮ‬ vs. ܶ. Any anomalies in the underlying TSL physics models as a function of ܶ could be easily identified, which is difficult or impractical with critical benchmarks. As light water is a very important neutron moderator for reactor and criticality safety applications, this physics-based validation approach provides a valuable tool to investigate and address any deficiencies in the TSL models and ensure the best quality of data is available to end users.