Radiation Transport in Random Media With Large Fluctuations

Neutral particle transport in media exhibiting large and complex material property spatial variation is modeled by representing cross sections as lognormal random functions of space and generated through a nonlinear memory-less transformation of a Gaussian process with covariance uniquely determined by the covariance of the cross section. A Karhunen-Loève decomposition of the Gaussian process is implemented to efficiently generate realizations of the random cross sections and Woodcock Monte Carlo used to transport particles on each realization and generate benchmark solutions for the mean and variance of the particle flux as well as probability densities of the particle reflectance and transmittance. A computationally efficient stochastic collocation method is implemented to directly compute the statistical moments such as the mean and variance, while a polynomial chaos expansion in conjunction with stochastic collocation provides a convenient surrogate model that also produces probability densities of output quantities of interest. Extensive numerical testing demonstrates that use of stochastic reduced-order modeling provides an accurate and cost-effective alternative to random sampling for particle transport in random media.


Introduction
Strongly heterogeneous media arise in several applications that include radiation shields, nuclear fuel, BWR moderators, clouds, planetary and stellar atmospheres, turbulent gases and plasmas.Neutral and charged particle transport computations in such media have relied heavily on formulating transport equations with spatially random coefficients (physical data) and developing solution methods to deal with the additional stochastic dimensions.Attempts at developing approximate closures that yield only low order statistical information (e.g., mean and variance of the flux) have proved to be highly restrictive under real physics conditions [1,2] or rely on techniques that require fluctuation amplitudes to be small for robustness [3].Recently, stochastic spectral methods such as polynomial chaos and stochastic collocation [4] have been developed for aleatoric uncertainty quantification and sensitivity analysis, and successfully applied in radiation transport work [5,6].Advances in UQ techniques have tended to focus on efficiently handling large numbers of uncertain variables but the rigorous stochastic basis of the approach also promotes its use in situations where stochasticity is due to spatial heterogeneity and the associated uncertainty is large.Here we apply these techniques to radiation transport in media with spatially randomly varying cross sections without restriction on fluctuation amplitudes.Specifically, we represent the cross sections as a lognormal spatial random process with specified mean, variance and covariance function and use a KL-decomposition to efficiently generate cross section realizations that are strictly positive.Woodcock Monte Carlo [7,8] is then used to simulate transport using random sampling of cross sections and deterministic sampling based on a stochastic collocation technique.A polynomial chaos expansion (PCE) model is built for the stochastic response, which enables construction of probability density functions of physical quantities of interest.Numerical results for the mean and variance in the scalar flux and leakage currents are obtained.We focus in this work on demonstrating our approach to solving random medium transport problems and defer application to specific problems to a future investigation.

Problem Statement and Stochastic Model
The random transport equation of interest is given by μ ∂ψ(x, μ, ω) ∂x where x and μ are the particle spatial and angular variables and the label ω denotes an element in the sample space, i.e., a particular realization.Macroscopic total or scattering cross section r, r = {t, s}, is represented as the product of microscopic cross section r and the stochastically dependent material atom density N at : The atomic density N at (x, ω) is modeled as a lognormal random field derived from Gaussian random field g(x, ω), yielding all macroscopic cross sections for the material and ensuring that these values are always positive.

The Karhunen-Loève Expansion
The optimal spectral expansion that maintains secondorder statistics of the Gaussian random field g(x, ω) is given by the Karhunen-Loève (KL) expansion, where g and v g are the average and variance of the Gaussian process and γ k , u k , and ξ k are the eigenvalue, eigenfunction, and random variables associated with the KL expansion.The KL eigenvalues and eigenfunctions are given by the solution of the following Fredholm integral equation of the second kind with the normalized covariance as the kernel [9]: where c g (x, x ) = C g (x, x )/v g , and C g (x, x ) is the process covariance.The eigenvalues are all positive and ordered with decreasing magnitude, which makes computation feasible by truncation of the KL expansion at an order that captures most of the variability, depending on the strength of correlation [5,9].The eigenspectrum can be obtained analytically for several covariance functions including an exponential covariance, characterized by a correlation length λ c , and linear covariance [9].The memory-less nonlinear transformation of the Gaussian process requires solving for the mean, variance, and covariance of the underlying Gaussian process in terms of the statistics of the cross section that is modeled as a lognormal process.It is not hard to relate the mean and variance of the processes, and the normalized covariances are related [5] by The eigenspectrum of the Gaussian KL process must be solved numerically since the Fredholm equation, Eq. ( 5), is not solvable analytically for the covariance of the Gaussian process given by Eq. ( 7).

Numerical Solution of the KL Eigenspectrum
The Nyström method numerically solves the Fredholm integral equation for the Karhunen-Loève eigenspectrum by discretizing over the domain using nodes x j and weights w j , where N Ny is the number of nodes in the discretization, and eigenvector values are solved at nodes x i .The numerical solve of the Fredholm equation is written in matrix notation as where matrix C is a symmetric positive semi-definite matrix with elements C i j = c g (x i , x j ), matrix W is a diagonal matrix with values w j , and vector u k spans the physical domain for each KL eigenvalue k.Eq. ( 9) is reformed as where the matrix B and vector u * k are given by: and W − 1 2 is the diagonal matrix with values √ w j .The eigevector can then be evaluated as: The simplest way to evaluate eigenfunction u k (x) using eigenvector u k is to return the eigenvector value which corresponds to the eigenvector node x i closest to x.One evaluation of u k (x) requires sampling only one value of u k and the resulting eigenfunction is a discontinuous step function in physical space.The accuracy of u k (x) is improved by evaluation using an interpolating function provided in the Nyström method theory [10,11]: The resulting eigenfunctions u k (x) are more accurate than when using the first interpolation scheme, continuous, and differentiable, though sampling a value of x requires N Ny evaluations of u k , making this interpolation scheme computationally expensive.A third interpolation scheme linearly interpolates between values of the eigenvector u k and defaults to use of the interpolating function in Eq. (13) when x is near the edge of the slab and therefore not bounded by values of u k .This interpolation scheme appears to be nearly as accurate as the scheme given in the Nyström theory, is continuous but not differentiable, and assuming that each value of x is equally likely to be sampled has a cost of 1 + 2 N Ny −1 N Ny -roughly three-evaluations of u k to sample a value of u k (x).This linear interpolation scheme is used in this work.
With the Gaussian process completely characterized and realizations of independent standard normal random variables readily sampled, realizations of cross sections are obtained from the truncated KL expression: To demonstrate the positivity and moment preserving properties of these lognormal processes, we plot fifty Gaussian random process realizations and fifty lognormal random process realizations each with g = 0.75 cm and v g = 0.225 cm 2 over a slab of length L = 5 cm in Fig. 1 and Fig. 2 respectively.Where realizations of the Gaussian  process contain negative values, realizations of the lognormal process do not.In order to recover the input variance, the lognormal process contains more values near zero and fewer values which may be significantly larger than the process mean.
As an implementation verification test and to demonstrate the increased accuracy of the Karhunen-Loève expansion with larger truncation orders K, the mean and standard deviation of an ensemble of 50,000 realizations were sampled at 20 different points in a slab using KL truncation orders of 3, 7, and 15.The observed sample mean and standard deviation values are plotted in Fig. 3 along with the input mean ( Σ r = 0.75 cm) and standard deviation ( shapes of the eigenfunctions can be seen in the oscillatory shape of the sampled moments, which are mitigated by increasing number of terms retained in the KL-expansion.

cm). Artifacts of the spatial
The moments of the cross sections are preserved over the ensemble with more accuracy for larger truncation orders K.

Monte Carlo Transport Solver for Spatially Continuous Cross Sections
Numerical simulation of particle transport on different realizations of the stochastic transport equation can be performed using a variety of methods.Deterministic solution methods require discretization of the geometry and, unless finite elements are used with mesh boundaries and basis functions built to match the Nyström interpolating scheme, introduce an error source.The traditional approach to Monte Carlo transport requires spatially discrete cross section values which could be achieved through a discretization in physical space, again introducing an error.Monte Carlo transport could be effected exactly by computing distance to collision based on integrals of the total cross section, the optical thickness, though the cross section profiles of interest here are not analytically integrable requiring numerical integration which introduces error and is likely to be expensive.Transport is solved here introducing only resolvable statistical error using Monte Carlo transport with Woodcock sampling.Monte Carlo transport with Woodcock sampling [7], Woodcock Monte Carlo (WMC) for short, samples a distance to potential collision based on a fictitious total cross section equal to or larger than the maximum total cross section over a particle flight.The value of the total cross section is compared to the value of the fictitious total cross section using a pseudo-random number to determine if the particle undergoes a collision or continues streaming.As long as WMC computations are not dominated by rejected potential interactions, WMC efficiently and exactly, within statistical noise, computes transport results.We note that WMC is used regularly in at least one production Monte Carlo transport code [7].The WMC solver used in this work has been utilized previously [8] and benchmarked against a suite of problems for which an analytic or near analytic solution is available.

Stochastic Solution Methods
The Karhunen-Loève expansion with a truncation of order K is used to characterize a random field as a function of K Gaussian random variables.Stochastic solution methods introduced here resolve the effects of these remaining random variables yielding statistical averages of particle flux and current and probability density function of leakage currents.
Statistical averages of a quantity of interest ϕ, here either particle flux in a cell or a leakage current, are defined with respect to the joint probability density P(ξ 1 , . . ., ξ N ) over all random variables retained in the KL-expansion as (15) Standard normal random variables in the KL expansion are known to be independent [9] such that the joint probability density factorizes into a product of standard normals for each random variable, vis., The integral in Eq. ( 15) is solved in one of two ways here, either through random or deterministic sampling of the random variables.Random sampling solutions provide probability density functions of the response built from a number of values equal to the number of transport evaluations.Solutions from the deterministic sampling approach are used to create a surrogate model of the response from which probability density functions are constructed using an arbitrary number of samples.Once a surrogate model is built, generation of probability density functions is much cheaper than generation using random sampling.

Random Sampling
Moments of the response can be solved by numerically evaluating the integral on the right-hand-side of Eq. (15) using random sampling (RS): where R is the total number of sampled realizations and ξ (i) k are the randomly sampled values of the standard normal random variables.Transport results are obtained for each of R realizations to yield ensemble averaged response moments which converge towards the solution as R −0.5 .The approach is dimensionally agnostic, meaning the convergence rate is not affected by the number of random dimensions.This method is appropriately called a Monte Carlo method, but we reserve that title in this work for Monte Carlo particle transport to avoid confusion.

Stochastic Collocation
Moments of the response can also be solved by numerically evaluating the integral on the left-hand-side of Eq. ( 15) using multi-dimensional quadrature: where Q k is the order of quadrature used on the k-th random variable and w (q k ) k and ξ (q k ) k are the weights and nodes of quadrature points.Such application of quadrature over the stochastic domain is called stochastic collocation (SC).Gauss-Hermite quadrature, being the optimal non-nested quadrature for integration of polynomials over Gaussian probability basis functions, is used, and if the response surface is sufficiently smooth, such Gauss-Hermite quadrature is effective.The stochastic collocation method, for isotropic collocation grids, Q k = Q k ∀k, k , convergences exponentially in log-log space.It does, however, suffer from the curse of dimensionality such that its convergence is slower with an increased number of random variables.The number of random variables below which SC is more efficient than RS depends on the smoothness of the response surface and the error tolerance desired in the response-smoother surfaces and tighter error tolerances favor SC over RS.For reasonably smooth problems and typical error tolerances, SC tends to be more efficient than RS for roughly 10 to 20 random variables.For the tensor product collocation grids used in this work the total number of realizations over which transport must be solved is equal to the product of the quadrature order in each dimension: Anisotropic stochastic collocation grids, for which ∃k, k : Q k Q k , are especially effective when integrating over KL random variables since KL terms are known to decrease monotonically in magnitude [12], though it is difficult to a priori know the optimal quadrature orders Q k .Adaptive methods which estimate the error contribution of each dimension may be useful in choosing the optimal order.Additional efficiency gains may be achievable using sparse collocation grids [5,6,12].

The Polynomial Chaos Expansion
While the stochastic collocation method does not provide full distributions of response values, the SC method can be used to build a surrogate model of the response using a polynomial chaos expansion (PCE).The polynomial chaos expansion using Hermite polynomials can be expressed as: where instances of u represent coefficients to be determined and H n represent n-order multivariate probabilists' Hermite polynomials.The K-dimensional PCE truncated at a total polynomial order of I can be written more succinctly using multi-indices i and ξ as where the total number of PCE terms is equal to Since the probability densities in each dimension are independent, the multidimensional Hermite polynomials factor into the product of univariate Hermite polynomials: The coefficients of the PCE are solved by projecting the expansion in Eq. ( 21) over the polynomial chaos and using stochastic collocation to compute the resulting inner products: where the angular brackets represent integration over the phase space with weight function given by the probability density p( ξ): The right-hand-side of Eq. ( 24) is simplified using orthogonality of the Hermite polynomial: to obtain for the expansion coefficients: where a i is defined for one dimension in Eq. ( 26).The integral on the right-hand-side of Eq. ( 27) is then numerically evaluated using the stochastic collocation weights w (q k ) k , nodes ξ (q k ) k , and solutions ϕ( ξ (q) ) generated when solving response moments in Eq. ( 18): H i 1 (ξ (q 1 ) 1 ) . . .
where q is a multi-index of values q k .The polynomial order in any stochastic dimension should not be larger than the quadrature order used in that dimension.In this work we set the total PCE polynomial order equal to the largest quadrature order used, vis., but additionally truncate any terms for which the polynomial order in a given dimension is greater than the quadrature order used in that dimension.The resulting number of PCE terms is thus The model constructed for a KL truncation order K will be more accurate for higher polynomial orders I and collocation orders Q k .Once constructed the surrogate model can be used to approximate the response corresponding to a new set of random variable values without solving transport over an additional realization.The surface can be cheaply sampled to, for instance, produce a visual of the response surface, estimate the response for a new set of input parameters, or create a probability density function of the response.
A polynomial fit of the system response does not guarantee purely positive values, especially when the response is near zero over much of stochastic domain.We introduce a positive-preserving transformation of the PCE by fitting the response to where n is a free parameter.Coefficients are then solved by applying stochastic collocation over the integral of A value of n = 0.5 has no effect on evaluation of coefficients for transport problems since transport response is always positive, and simply takes the absolute value of PCE values.Values of n smaller than 0.5 map the response values less than 1 to smaller positive values and values greater than 1 to larger values, creating a more sharply varying surface, whereas values of n greater than 0.5 map all values closer to 1.We hypothesize that mapping values which the PCE will fit closer to 1 increases the accuracy of the fit, as long as values of n are not chosen to be so large that numerical roundoff of exponential computations introduces significant error.As the value of n approaches infinity, all values map to 1 within machine precision and all information defining the surface is lost.We somewhat arbitrarily have chosen a value of n = 1 for the results in this paper.Other positive-preserving transformations are certainly conceivable and are a topic of future work.
With or without a positive-preserving transformation, finite polynomial fits of a response which spans an unbounded domain will, in regions of low probability, include values which are very high.The simplest method for handling such very large sampled values when constructing a probability density function is to ignore them since they are sampled infrequently.Perhaps a screening method which rejects samples taken from regions of sufficiently low probability would treat this problem.

Numerical Results
Solution of the stochastic transport equation is demonstrated on a problem described by microscopic total and scattering cross sections σ t = 2.0 cm 2 and σ s = 1.0 cm 2 , atomic density average and variance of N at = 1.0 cm −3 and v N at = 1.0 cm −6 , and an exponential covariance function with correlation length λ c = 1.5 cm.
Realizations of the random cross sections are generated using a lognormal transformation of the Karhunen-Loéve (KL) expansion truncated at expansion order K = 7 in which the eigenspectrum is solved using the Nystöm method (NM) with an N Ny = 100 node discretization.Transport results are obtained using Woodcock Monte Carlo (WMC) with N = 1.6 × 10 5 particle histories.KL random variable values are chosen either using random sampling (RS) with R = 2 × 10 5 samples, isotropic stochastic collocation (iSC) with quadrature orders Q = {5, 5, 5, 5, 5, 5, 5} for a total of R S C = 78, 125 realizations, or anisotropic stochastic collocation (aSC) with quadrature orders Q = {6, 5, 5, 4, 4, 3, 3} for a total of R S C = 21, 600 realizations.Probability density functions (PDFs) of the reflectance and transmittance are generated using R = 2 × 10 5 samples when using the RS method, and by sampling from a constructed polynomial chaos expansion (PCE) model, both with and without the positivepreserving transformation, 1 × 10 6 times when solving using SC.
These solver parameters were chosen as the result of convergence studies and the transmittance standard deviation relative error is estimated using these parameters to be less than 3 %.Perhaps the most noteworthy convergence study was the one performed on the choice of Karhunen-Loève truncation order K, using the transmittance standard deviation as the output of interest.A KL truncation order of K = 10 was used to estimate the correct solution, holding all other solver parameters the same, and the relative difference of the solution given by values of K ≤ 10 and K = 10 is plotted in Figure 4.The estimated relative error contribution to the transmittance standard deviation based on the truncation order of the KL expansion is therefore 0.8 %.A similar study was performed for choosing  the number of particle histories N, the number of random samples R when using random sampling, and the quadrature orders Q when using stochastic collocation.The SC convergence study was performed using isotropic stochastic collocation such that the quadrature order in each dimension was increased at that same rate.Error was estimated using several anisotropic SC grids chosen based on intuition, and one of them was chosen for this problem.Problem parameters were chosen for each solver such that the estimated bound on the total error of the transmittance standard deviation was about 3 %.
Table 1 records the transmittance and reflectance values yielded by the RS, iSC, and aSC methods and the number of realizations used to produce these results.The greatest difference in transmittance standard deviation values between these solutions is about 1.6 %, less than the estimated bound on the error of 3 %.The largest relative difference between results for any quantity is a relative difference of about 5.6 % for the reflectance standard deviation.Error in these final quantities can be reduced by refining solver parameters such as choosing a larger KL truncation order K, using more particle histories, or increasing the number of realizations either directly in random sampling or by using higher-order quadrature in the SC method.Data from convergence studies, such as the plot in Figure 4, can be used to estimate the change of system parameters that will reduce the total error for the lowest computational cost.
Having solved the problem using each method to roughly the same accuracy, it is clear that for this problem the anisotropic stochastic collocation implementation solves the problem most efficiently, requiring about 10 times less realizations than random sampling.There is likely a more efficient set of anisotropic quadrature orders than the somewhat arbitrarily chosen set used here and an appropriate adaptive method may be able to find that set of quadrature orders efficiently.Isotropic stochastic collocation is about 2.5 times more efficient than random sampling for this problem.The SC methods are expected to be yet more efficient for more highly correlated problems in which the correlation length λ c is larger, since the magnitude of successive KL terms falls off more quickly for more correlated processes allowing more liberal truncation of the expansion and a higher degree of anisotropy.
Scalar flux mean and relative standard deviation solved in 100 cells throughout the slab are given in Figures 5  and 6.Scalar flux values yielded by each of the three methods match up well with one another.The greatest deviation between the methods appears to be near the transmissive boundary at which the mean flux is a little higher and relative standard deviation is a little lower for the isotropic stochastic collocation method than the others.This observation is consistent with the transmittance values reported in Table 1.
Whereas the mean and standard deviation of leakage values provide information often useful in shielding applications, information beyond the first two moments may be of value.For example, the median expected transmittance or the largest credible transmittance may be of more interest.Probability density functions are therefore generated which provide a variety of information beyond that given by the mean and standard deviation of leakage values.
Probability density functions of the transmittance generated with random sampling and isotropic stochastic collocation are given in Figure 7 and with random sampling and anisotropic stochastic collocation in Figure 8.Each plot contains results generated by constructing a PCE model in a straightforward way and using the positivepreserving transformation introduced here with free parameter n = 1.PDFs generated using the transformation with a free parameter value of n = 0.5 would be the same as those generated by the PCE model without transformation but in which all negative values are tallied as positive values of the same magnitude.Five times as many sam- ples were taken to form the PDFs constructed from PCE models than those constructed using RS since it is much cheaper to sample a PCE model of the response once constructed than to solve transport results over additional realizations using RS.The PDF generated from the naive implementation of the PCE contains negative values due to the truncation of the PCE, making the positive-preserving transformation useful to maintain physically meaningful results.With or without this transformation, the polynomial fit over an unbounded support given by the PCE contains a low probability of sampling artificially large values.A method to screen out such samples may be worth investigating.The tapering off of high values in the PDF produced using iSC is also likely due to use of a polynomial fit over an unbounded support.It appears that having a quadrature order of 6, as in the anisotropic case, as opposed to 5, as in the isotropic case, in the first KL term re- duces much of this effect.The transmittance PDFs generated from PCE models using the PCE transformation presented here contain less larger values.Comparing against the PDF produced using RS, this appears to be a good trait in the iSC case for this problem and a bad trait in the aSC case for this problem.
The PDFs produced using the PCE model closely match the PDFs generated with random sampling, and provide more than a mean and It is unclear whether using the positive-preserving transformation produced PDFs that were more accurate than the PDFs produced without the transformation on the positive portion of the domain, but the positive-preserving transformation did a good job a resolving the steep drop-off in transmittance values near zero.It is possible that other positive-preserving transformations, such as a maximum entropy representation, may prove to be a better choice than the one presented here; this is a topic of future work.
Likewise, PDFs of the reflectance generated using random sampling and isotropic stochastic collocation are given in Figure 9 and PDFs of the reflectance generated using random sampling and anisotropic stochastic collocation are given in Figure 10.While maintaining posi-  tive values did not turn out to be as important in modeling the reflectance response as modeling the transmittance response, the results generated using the positive-preserving transformation do not appear less accurate than when not using the transformation.PDFs of the response generated from PCE models deviate somewhat more from the PDF generated using random sampling for the reflectance than the transmittance.This is not an unexpected result for the problem parameters chosen since the standard deviation of the transmittance was the quantity that solver parameters were chosen to reduce the error of below a certain tolerance and the relative difference in the reflectance values were larger than the relative difference in transmittance values (Table 1).The PDFs generated from the PCE models could be made more accurate by increasing one or more of the solver parameters such as using more KL terms, more particle histories, more random samples, or a higher PCE truncation yielding more higher-order polynomial cross terms in the expansion.

Conclusions
We have demonstrated that a combination of stochastic spectral representation and Woodcock Monte Carlo simulation allows efficient computation of radiation transport in strongly random media.The use of a lognormal distribution to represent the fluctuations in medium properties facilitates enforcement of strict positivity of cross section realizations, a condition that has challenged standard approaches that assume Gaussian fluctuations [8].The fact that the lognormal random process is a memory-less nonlinear transformation of a Gaussian process further enables efficient and accurate reconstruction of the random cross section realizations using a Karhunen-Loève representation of the cross section.The robustness of this model enables fluctuations of arbitrarily large amplitude to be studied, using both random sampling of Karhunen-Loève realizations and more efficient stochastic collocation techniques.The polynomial chaos expansion is shown to be a useful surrogate model of the response populated using stochastic collocation results which allows for efficient construction of probability density functions of out-put variables.The use of a positive-preserving transformation of the polynomial chaos expansion further increases its usability by enforcing that all output variable values are positive with no apparent loss of accuracy in the representation.Further work on these techniques may include examination of covariance functions other than the exponential covariance, advanced stochastic collocation grids such as adaptive or sparse grids, use of a different positivepreserving transformation for the surrogate model of the response such as a lognormal transformation, or more judiciously selected ceiling cross sections to enhance the efficiency of the Woodcock Monte Carlo solver.

Figure 3 .
Figure 3. Observed Mean and Standard Deviation in Ensemble of Lognormal Random Cross Section

Figure 4 .
Figure 4. Standard Deviation Relative Error of KL Truncation

Figure 5 .
Figure 5. Mean Scalar Flux in Slab

Table 1 .
Leakage Values and Cost