Experimental and Simulated Spectral Gamma-Ray Response of a NaI(Tl) Scintillation Detector used in Airborne Gamma-Ray Spectrometry

The objective of this work is to simulate the spectral gamma-ray response of NaI(Tl) scintillation detectors for airborne gamma-ray spectrometry (AGRS) using the state-of-the-art multi-purpose Monte Carlo code FLUKA. The study is based on a commercial airborne gamma-ray spectrometry detector system with four individual NaI(Tl) scintillation crystals and a total volume of 16.8 L.

To validate the developed model, radiation measurements were conducted using ⁵⁷Co, ⁶⁰Co, ⁸⁸Y, ¹⁰⁹Cd, ¹³³Ba, ¹³⁷Cs and ¹⁵²Eu calibration point sources with known activities and source-detector geometries under laboratory conditions. In addition, empirical calibration and resolution functions were derived from these measurements combined with additional radiation measurements adopting natural uranium, thorium and potassium volume sources.

The simulation results show superior accuracy and precision compared to previous AGRS simulation models with a median relative spectral error < 10 % for most of the radiation sources. Moreover, the implementation of a lower level discriminator model and detailed modelling of the laboratory result in a significant improvement in model accuracy for spectral energies < 100 keV compared to previous studies. Yet thorough statistical analysis incorporating statistical and systematic uncertainty estimates revealed statistically significant deviations between the simulated and measured spectra in the spectral region around the Compton edge, which could be attributed to the scintillator non-proportionality.

These findings imply that the linear energy deposition model applied in this and previously developed AGRS simulation models should be revised and considered to be replaced by more accurate non-proportional models.

1 Introduction

Airborne gamma-ray spectrometry (AGRS) is an established technique to identify and quantify terrestrial radionuclides using a gamma-ray spectrometer mounted in an aerial vehicle (Erdi-Krausz et al., 2003). This technique was originally developed in the early second half of the past century for radioactive mineral exploration and geological mapping (Grasty, 1975). In the years that followed the first successful geological and mineral exploration surveys, AGRS has proven to be also an essential tool in a variety of other fields, including emergency response to severe nuclear accidents, orphan source location, non-radioactive mineral exploration, nuclear facility monitoring, soil mapping, snow coverage quantification, geothermal exploration and landslide mapping (Connor et al., 2016; Li et al., 2020; Pradeep Kumar et al., 2020). The severe nuclear accident at the Fukushima Daiichi nuclear power plant in 2011 (Torii et al., 2013) as well as the concerns regarding nuclear terrorism prompted new research efforts leading to the development of new innovative AGRS systems, especially in the field of unmanned aerial vehicles (UAV) and directional-sensitive-based AGRS systems (Connor et al., 2016; Kulisek et al., 2018).

As pointed out by Fortin et al. (2017), despite decades of development, there are still some unsolved problems in the processing as well as in the evaluation of AGRS data. Due to the complexity of the radiation transport and the detector response, the calibration of AGRS systems, i. e. the derivation of a reference system to quantify the terrestrial radionuclide concentrations based on the spectral detector response, remains a challenge.

Up to now, experimental calibration procedures based on calibration pads are the suggested approach according to the technical guidelines published by the International Atomic Energy Agency (IAEA) (Erdi-Krausz et al., 2003). However, these calibration-pads-based methods have two main limitations. First, the required calibration pads are expensive to build as well as to decommission. Second, due to the required source activities and radiation protection legislation, only a very limited range of radionuclides is available for experimental calibration purposes of AGRS systems.

With increasing computational capabilities and development of high-fidelity Monte Carlo based radiation transport codes over the past two decades, simulation-based calibration approaches are now commonly used for laboratory and in-situ spectrometry (Van der Graaf et al., 2011), in the marine environment (Androulakaki et al., 2016) or in space exploration (Prettyman et al., 2006). Despite the apparent advantages, i. e. unrestricted capabilities in terms of source characteristics and detector modelling as well as cost effective and fast evaluation, simulation-based calibration has not established itself as an equivalent methodology for standard AGRS systems, which are defined in this study as manned aircraft systems using a spectrometer with multiple 10.2 cm × 10.2 cm × 40.6 cm prismatic NaI(Tl) scintillation crystals.

First attempts to simulate the detector response of a standard AGRS system were performed by Allyson and Sanderson (1998) using an in-house Monte Carlo code. Significant systematic deviations (> 10 %) between measurement and simulation results were found and attributed to the simplified model, i. e. the multi-component detector system was simplified to a single crystal together with an equivalent aluminum casing. Billings and Hovgaard (1999) adopted a modified version of the Monte Carlo code developed by Allyson and Sanderson (1998) and performed extensive simulations to characterize the angular detector response of a standard AGRS detector system. Again, the multi-component detector system was simplified to a single crystal. In addition, only the scintillation crystal was incorporated in the model, i. e. all casing materials were neglected, and no spectral analysis was performed. Torii et al. (2013) adopted the EGS5 Monte Carlo code to calibrate their standard AGRS system for the aerial radiological surveys, which were performed to quantify the dispersion of ¹³¹I and ¹³⁴Cs in the aftermath of the nuclear accident in the Fukushima Daiichi nuclear power plant. Similar to Billings and Hovgaard (1999), Torii et al. (2013) did not perform spectral analysis but used the simulation to derive some empirical parameters for their analytical mono-energetic model. Sinclair et al. (2016) utilized the EGSnrc Monte Carlo code to calibrate their standard AGRS system for the quantification of the radionuclide ¹⁴⁰La. Again, only peak count rate measures were extracted from the performed simulations. In contrast to the earlier studies, Torii et al. (2013) as well as Sinclair et al. (2016) kept the multi-component crystals and included also additional detector components such as insulation materials, photomultiplier tubes (PMT) and multi-channel analyzers (MCA) in their simulation models. Zhang et al. (2018) adopted the MCNP Monte Carlo code to simulate the spectral detector response of a standard AGRS system. They included primitive cubic shells to model the aircraft but did not provide any evidence for model validation using measurement data.

As it emerges from this introduction, studies covering simulation-based calibration of standard AGRS systems are scarce. Moreover, most of the studies focus on full energy peak characterization for single radionuclides, adopted significant geometrical as well as physical simplifications, did not assess systematic uncertainties and, for some cases, lack experimental validation.

To address this gap, we performed high-resolution Monte Carlo simulations of a standard AGRS system under laboratory conditions using the state-of-the-art radiation transport code FLUKA (Böhlen et al., 2014; Battistoni et al., 2015). The detector model includes all relevant components such as scintillation crystals, casings, insulation materials, PMT and MCA in high detail. In addition, benchmark measurements were performed with various radionuclides adopting a centered source-detector configuration in the calibration laboratory at the Paul Scherrer Institute. As a first step towards a simulation-based calibration methodology, this study is concerned with centered source-detector configurations and not with angular detector response analysis. By comparing the experimental and simulated spectral results, new insights in the current understanding of the detector response can be obtained. Furthermore, adopting rigorous uncertainty quantification, the statistical as well as systematic uncertainties of the current model can be assessed and included in a statistically sound interpretation of the simulation results. For the reader's convenience, a full list of all symbols and acronyms used in this study is attached in Appendix C.

2 Methods

2.1 Radiation measurements

The goal of the laboratory-based radiation measurements are twofold. First, the spectral measurements are used to derive crystal-specific energy calibration and detector resolution parameters. Second, the radiation measurements can be adopted to validate the simulated spectral detector response and thereby quantify the model accuracy and precision.

2.1.1 Experimental setup

This study focuses on a commercial AGRS system, which is used regularly in annual survey campaigns by civil and military institutions in Switzerland for nuclear facility monitoring, geophysical studies and radiological emergency response training (Butterweck et al., 2021). This standard AGRS system was developed by Mirion Technologies Inc., France, and consists of four 10.2 cm × 10.2 cm × 40.6 cm prismatic NaI(Tl) scintillation crystals protected by individual aluminum casings together with related PMT and MCA. The scintillation crystals as well as the associated electronics are embedded in a thermal-insulating and vibration-damping polyethylene (PE) foam and installed in an aluminum box as illustrated in Fig. 1c. The spectrometer features 1024 channels for an energy range between 30 and 3000 keV together with automatic linearization of the individual scintillation crystal spectra.

Figure 1 Experimental setup in the calibration laboratory at the Paul Scherrer Institute (PSI). (a) Radiation source including the source holder, the radiation source disk (25 mm × 3 mm) and the ion exchange resin sphere (1 mm), in which the activity mass is embedded. (b) Cut view of the laboratory room (5.3 m × 4.5 m × 3 m) with the detector mounted on an aluminum frame and the source on a tripod. (c) Cut view of the detector box (86 cm × 60 cm × 30 cm) including the internals such as the scintillation crystals, photomultiplier tubes (PMT), multi-channel analyzers (MCA) and polyethylene (PE) foam. All figures were created using the graphical interface FLAIR (Vlachoudis, 2009). For better visibility and interpretability, false colors are adopted.

Radiation measurements where performed in the calibration laboratory at the Paul Scherrer Institute in Switzerland using standard calibration disk sources from the Eckert & Ziegler Nuclitec GmbH (⁵⁷Co, ⁶⁰Co, ⁸⁸Y, ¹⁰⁹Cd, ¹³³Ba, ¹³⁷Cs and ¹⁵²Eu). These sources consist of a 3 mm thick plastic disk with a 25 mm diameter, in which the active source mass is embedded in a 1 mm diameter ion exchange resin sphere (cf. Fig. 1a). Specific source properties including the activity at the corresponding measurement date, half life and major photon emission lines together with corresponding uncertainty values are summarized in Table A1 in Appendix A. Additional volume sources with unknown activities were used for energy and spectral resolution calibration, i. e. natural potassium K_nat, natural uranium U_nat and natural thorium Th_nat (cf. Sect. 2.1.2).

The detector box was placed on an aluminum frame and vertically oriented along the detector axis of symmetry, i. e. the z-axis, for all laboratory-based radiation measurements as illustrated in Fig. 1b. The standard calibration disk sources were inserted in a low absorption source holder and mounted on a tripod (cf. Fig. 1a). The source holder was specifically developed and constructed for these measurements using additive manufacturing techniques. The ion exchange resin sphere inside the calibration disk source was then aligned on the detector x-axis and placed in (1.00 ± 0.01) m distance to the detector box, where the outer surface of the detector box was defined as the reference surface. The detector x-axis is centered on the scintillation crystals as shown in Fig. 1c. To measure the source-detector distances and to position the sources accurately, distance as well as positioning laser systems were used. Between radiation measurements, background measurements were performed regularly for background correction and gain stability checks. For all measurements, the air temperature as well as the air humidity in the calibration laboratory was controlled by an air conditioning unit and logged by an external sensor. The air temperature was set at (18.8 ± 0.4) ^∘C and the relative air humidity at (42 ± 3) %. The ambient air pressure, which was also logged by the external sensor, fluctuated around (982 ± 5) hPa. If not otherwise stated, uncertainties are provided as 1 standard deviation (SD) values in this study (coverage factor k = 1).

During measurements, additional instruments and laboratory equipment were located in the calibration laboratory, e. g. shelves, a workbench, a source scanner or a boiler as shown in Fig. 1b. The effect of these features on the detector response is carefully assessed in a sensitivity analysis discussed in Sect. 2.2.1.

2.1.2 Spectral postprocessing

As a first step in the postprocessing of the measured pulse-height spectra, the gross count spectra C_gr,k are transformed to background and dead-time-corrected net count rate spectra c_exp,k normalized by the source activity A by subtracting the background count rate spectrum c_bg,k from the gross count rate spectrum c_gr,k for each channel k:

$$\begin{array}{}\text{(1)}& c_{\exp ,k}={\displaystyle \frac{\mathrm{1}}{A}}\left({\displaystyle \frac{C_{\mathrm{gr},k}}{t_{\mathrm{gr}}}}-{\displaystyle \frac{C_{\mathrm{bg},k}}{t_{\mathrm{bg}}}}\right)=c_{\mathrm{gr},k}-c_{\mathrm{bg},k}\end{array}$$

The dead time corrected background and gross count rate spectra c_bg,k and c_gr,k are obtained by dividing the corresponding count spectra C_bg,k and C_gr,k by the recorded measurement live time t_bg and t_gr, respectively. The resulting net spectra are computed in units of counts per second [cps] and per source activity [Bq], i. e. [cps Bq⁻¹]. Additional information on the uncertainty analysis for the net count rate spectra c_exp,k as well as the source activity A can be found in Appendix A.

To characterize the relationship between the spectral energy and the detector channel as well as to quantify the detector resolution, statistical measures for the center and dispersion of the full energy peaks (FEP) in the pulse height spectra are required. In this study, the center and dispersion of the FEP are quantified by the mean μ and standard deviation σ of a Gaussian model (Knoll, 2010). The net count rate spectrum $c_{\exp ,k^{\prime }}$ containing a general multiplet of N_pk FEPs is then described as the sum of N_pk independent Gaussian peaks together with a baseline model $c_{\mathrm{b},k^{\prime }}$ over a predefined connected subset of channels $\mathit{\{}{k}^{\prime }:H_{k^{\prime }}\in \mathit{\{}\mathrm{1},\mathrm{\dots },\mathrm{1024}\mathit{\left\}}\mathit{\right\}}$:

$$\begin{array}{}\text{(2)}& c_{\exp ,k^{\prime }}=c_{\mathrm{b},k^{\prime }}+\sum _{l=\mathrm{1}}^{N_{\mathrm{pk}}}{a}_l\cdot e^{-\frac{{\left(H_{k^{\prime }}-{\mathit{\mu }}_l\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_l^{\mathrm{2}}}}\end{array}$$

where a_l is an additional normalizing model parameter for each FEP. For the baseline, a physics-oriented empirical model was adopted (Westmeier, 1992). Weighted non-linear least-squares (WNLLS) regression using the interior-reflective Newton method (Coleman and Li, 1996) was utilized to determine the 3N_pk statistical model parameters, i. e. a_l, μ_l and σ_l as well as corresponding standard deviation values ${\mathit{\sigma }}_{a_l}$, ${\mathit{\sigma }}_{{\mathit{\mu }}_l}$ and ${\mathit{\sigma }}_{{\mathit{\sigma }}_l}$ with $l\in \mathit{\{}\mathrm{1},\mathrm{\dots },N_{\mathrm{pk}}\mathit{\}}$. For the weights, the squared inverse of the coefficient of variation, i. e. $c_{\exp ,k^{\prime }}^{\mathrm{2}}/{\mathit{\sigma }}_{\mathrm{pois},\exp ,k^{\prime }}^{\mathrm{2}}$, was used (cf. Appendix A). In this study, the following emission lines were selected for FEP characterization using WNLLS: 1173.228(3), 1332.492(4) keV [⁶⁰Co]; 898.042(3), 1836.06(1) keV [⁸⁸Y]; 88.0336(1) keV [¹⁰⁹Cd]; 661.657(3) keV [¹³⁷Cs]; 121.7817(3), 344.279(1) keV [¹⁵²Eu]; 1460.822(6) keV [K_nat]; 2614.51(1) keV [Th_nat]; 351.932(2), 609.320(5) keV [U_nat], where 1 SD values are provided as least significant figures. Photon energies of those emission lines were taken from the ENDF/B-VIII.0 nuclear data file library (Brown et al., 2018). A representative singlet containing the FEP at 121.7817(3) keV originating from the ¹⁵²Eu source is shown in Fig. 2a illustrating the experimental net count rate spectrum together with the resulting baseline model, the WNLLS model as well as the 99 % prediction interval (PI). It is important to note that, in order to determine the Gaussian mean μ_l and standard deviation σ_l parameters of a multiplet, no activity normalization of the net count rate spectrum according to Eq. (1) is required.

Figure 2 Spectral postprocessing of the experimental data for the scintillation crystal 3. (a) Experimental net spectrum of the ¹⁵²Eu source showing the full energy peak at 121.7817(3) keV together with the resulting peak fit and baseline model. (b) Energy calibration model using peak fit data from various radiation sources, i. e. ⁶⁰Co, ⁸⁸Y, ¹⁰⁹Cd, ¹³⁷Cs, ¹⁵²Eu, K_nat, Th_nat and U_nat. (c) Experimental resolution model characterizing the full width at half maximum (FWHM) and adopting the same peak fit data from radiation sources as for the energy calibration model. The prediction interval (PI) is computed for a significance level of 1 % for all statistical models. Experimental error bars are defined as 1 SD (only visible in Fig. 2a).

2.1.3 Spectral energy calibration and detector resolution model

Using the Gaussian mean parameters μ_l from the previous spectral postprocessing step and the known true photon emission lines, a calibration regression model can be derived relating the spectral energy E with the detector channel H. Due to the spectral linearization performed by the spectrometer, a simple polynomial model of degree 1:

$$\begin{array}{}\text{(3)}& H=d_{\mathrm{1}}\cdot E\end{array}$$

with the slope parameter d₁ proved to be sufficient to describe the relation between the spectral energy E and the detector channel H adequately.

The quantification of the spectral resolution was achieved by deriving an analytical resolution model adopting an exponential function:

$$\begin{array}{}\text{(4)}& \mathrm{FWHM}\left(H\right)=b_{\mathrm{1}}\cdot H^{b_{\mathrm{2}}}\end{array}$$

where the full width at half maximum (FWHM) parameter is adopted in agreement with the standards in gamma-ray spectrometry (Knoll, 2010). The FWHM is analytically related to the Gaussian standard deviation σ_res by $\mathrm{FWHM}=\mathrm{2}\cdot \sqrt{\mathrm{2}\cdot \log \left(\mathrm{2}\right)}\cdot {\mathit{\sigma }}_{\mathrm{res}}$. In general, the resolution of scintillation detectors, specifically Na(Tl) crystals, can be attributed to a variety of factors, including quantum fluctuation in the number of information carriers, scintillator non-proportionality, crystal impurities or electronic noise (Knoll, 2010). Due to this complexity, there is no consensus in the literature on the functional form of an appropriate resolution model (Casanovas et al., 2012). The selection of the exponential model may seem therefore to a certain degree arbitrary but, due to its simplicity and versatility, has proven to be appropriate not only for this but also for other studies (Shi et al., 2002; Gardner and Sood, 2004).

For both models, WNLLS regression using the interior-reflective Newton method (Coleman and Li, 1996) was adopted with the squared inverse of the coefficient of variation, i. e. ${\mathit{\mu }}_l^{\mathrm{2}}/{\mathit{\sigma }}_{{\mathit{\mu }}_l}^{\mathrm{2}}$ and ${\mathit{\sigma }}_l^{\mathrm{2}}/{\mathit{\sigma }}_{{\mathit{\sigma }}_l}^{\mathrm{2}}$, as weights to determine the model parameters d₁, b₁, b₂ and corresponding standard deviation values ${\mathit{\sigma }}_{d_{\mathrm{1}}}$, ${\mathit{\sigma }}_{b_{\mathrm{1}}}$, ${\mathit{\sigma }}_{b_{\mathrm{2}}}$ for each scintillation crystal. Evaluated goodness-of-fit measures showed acceptable prediction properties, i. e. a coefficient of determination R² > 0.9999 and R² > 0.997 together with a relative generalization error < 0.02 % and < 0.2 %, for both the calibration and resolution model, respectively. To estimate the relative generalization error, the predicted residual error sum of squares, sometimes also referred to as the leave-one-out cross-validation error, normalized by the variance of the response variable was computed (Tarpey, 2000). In Fig. 2b–c, the resulting calibration and resolution model together with the experimental Gaussian mean μ_l and standard deviation σ_l parameters as well as the 99 % PI are shown for the scintillation crystal 3.

2.2 Detector response simulations

In this study, we adopt a single-stage stochastic model using the state-of-the-art multi-purpose Monte Carlo code FLUKA (Böhlen et al., 2014; Battistoni et al., 2015) version 4.2.1 together with the graphical interface FLAIR (Vlachoudis, 2009) version 3.1–15.1 to simulate the radiation transport as well as the energy deposition in the scintillation crystals. All simulations were performed on a computer cluster at the Paul Scherrer Institute utilizing parallel and high-performance computing.

2.2.1 Simulation model

In this study, the most accurate physics mode available in FLUKA, i. e. “precisio”, was used featuring a fully coupled photon, electron and positron radiation transport for our source-detector configuration. Lower transport thresholds were set to 1 keV for the scintillation crystals as well as the closest objects to the crystals, e. g. reflector, optical window and aluminum casing for the crystals. For the remaining model parts, the transport threshold was set to 10 keV to decrease the computational load while maintaining the high-fidelity transport simulation in the scintillation crystals. In accordance with previous studies (Allyson and Sanderson, 1998; Billings and Hovgaard, 1999; Sinclair et al., 2016; Kulisek et al., 2018; Zhang et al., 2018), a linear energy deposition model is applied in this study neglecting the scintillator non-proportionality of NaI(Tl) (Payne et al., 2009; Moses et al., 2012). For the radioactive sources, detailed photon, electron and positron emissions are simulated using the “raddecay” card in a semi-analogue mode according to the corresponding decay schemes allowing an event-by-event analysis.

It is important to add that for the ¹⁰⁹Cd source, a software bug was detected in the current FLUKA version and confirmed by the developers. Therefore, a custom source user routine was written using the emission data from the ENDF/B-VIII.0 nuclear data file library (Brown et al., 2018) to simulate the ¹⁰⁹Cd source adequately.

The energy deposition events in the scintillation crystals are scored individually on an event-by-event basis using the custom user routine “usreou” together with the “detect” card. The number of primaries N_pr was set to 10⁷ for all simulations, which guaranteed a maximum relative statistical standard deviation ${\mathit{\sigma }}_{\mathrm{stat},\mathrm{sim},k}/c_{\mathrm{sim},k}<\mathrm{1}\mathit{\%}$ and a maximum relative variance of the sample variance VOV_k < 0.01 % for all detector channels with ${\mathit{\mu }}_{\mathrm{LLD}}+\mathrm{2}\cdot {\mathit{\sigma }}_{\mathrm{LLD}}\le H_k\le d_{\mathrm{1}}\cdot E_{\mathrm{FEP},\max }+\mathrm{2}\cdot {\mathit{\sigma }}_{\mathrm{res}}$ and with respect to the summed detector crystal spectra (cf. Appendix A and B). For E_FEP,max, the maximum photon energy of the corresponding source using a detection threshold of 5 % with respect to the emission probability of the strongest line within the associated decay branch was used (cf. Table A1 in Appendix A).

The simulation model includes all relevant detector components such as scintillation crystals, reflectors, aluminum casings, insulation materials, PMT and MCA in high detail (cf. Fig. 1c). Geometrical as well as material properties were provided by the manufacturer of the Swiss AGRS system. The calibration disk sources as well as the source holder are modelled in high detail with known material and geometrical properties (cf. Fig. 1a). On the other hand, the laboratory room together with additional instruments and equipment, e. g. the workbench, shelves or a source scanner shown in Fig. 1b, are modelled in less detail. For this simplifications, care was taken to preserve the overall opacity as well as the mass density.

To assess the effect of the individual model elements on the detector response, a sensitivity analysis was performed changing one set of model elements at the time (OAT) for the ⁶⁰Co Monte Carlo simulation as a representative model. More specifically, starting with the high-fidelity model illustrated in Fig. 1 as the reference model (model A), specific element sets were excluded from this model and replaced with air. These element sets are defined as follows: laboratory room related elements {floor, walls, windows, door, ceiling} (model B), detector box related elements {aluminum box, aluminum frame, MCA, PMT, PE foam} (model C), source related elements {source disk, source holder, tripod} (model D) and additional laboratory equipment elements {shelves, workbench, source scanner, boiler} (model E).

In addition, there are some uncertainties regarding the elemental composition and mass density of the reflector, which is located between the aluminum casing and the scintillation crystals. Based on publicly available sources for our detector, the reflector is modelled with a polytetrafluoroethylene (PTFE) foil with an equivalent mass density of 2.25 g cm⁻³. In contrast, most of the other studies adopted MgO for the reflector with an equivalent mass density of 2 g cm⁻³ (Shi et al., 2002; Salgado et al., 2012) based on an early study by Saito and Moriuchi (1981). Consequently, an additional model (model F) was defined using MgO instead of PTFE for the reflector to quantify the sensitivity of this element on the detector response.

2.2.2 Spectral postprocessing

The obtained event-by-event energy deposition data is transformed to channel-based count rate spectra in three distinct steps for each of the four scintillation crystals. First, the individual energy deposition events are transformed from the spectral energy space E to the detector channel space H adopting the energy calibration models (Eq. 3) discussed in Sect. 2.1.3.

Second, to account for the finite detector resolution, the individual energy deposition events are broadened according to the derived resolution models for the individual crystals (Eq. 4). In more detail, the energy deposition events E_dep are assumed to be normally distributed with a mean and standard deviation equal to the simulated energy deposition E_sim and the corresponding standard deviation predicted by the resolution model σ_res, respectively:

$$\begin{array}{}\text{(5)}& E_{\mathrm{dep}}\sim \mathcal{N}\left(E_{\mathrm{sim}},{\mathit{\sigma }}_{\mathrm{res}}^{\mathrm{2}}\left(d_{\mathrm{1}}\cdot E_{\mathrm{sim}}\right)/d_{\mathrm{1}}^{\mathrm{2}}\right)\end{array}$$

The broadening, sometimes also referred to as Gaussian broadening (Van der Graaf et al., 2011), is performed in parallel with the binning and normalization of the recorded energy deposition events E_sim,i with $i\in \mathit{\{}\mathrm{1},\mathrm{\dots },N_{\mathrm{dep}}\mathit{\}}$ to generate the simulated count rate spectrum $c_{\mathrm{sim},k}^{*}$ for the detector channels H_k with $k\in \mathit{\{}\mathrm{1},\mathrm{\dots },\mathrm{1024}\mathit{\}}$ and the half channel width ΔH as follows:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{c}_{\mathrm{sim},k}^{*}& ={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{pr}}}}\cdot \sum _{i=\mathrm{1}}^{N_{\mathrm{dep}}}\underset{H_k-\mathrm{\Delta }H}{\overset{H_k+\mathrm{\Delta }H}{\int }}\mathit{\phi }\left({\displaystyle \frac{H-d_{\mathrm{1}}\cdot E_{\mathrm{sim},i}}{{\mathit{\sigma }}_{\mathrm{res}}\left(d_{\mathrm{1}}\cdot E_{\mathrm{sim},i}\right)}}\right)\mathrm{d}H\\ & ={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{pr}}}}\cdot \sum _{i=\mathrm{1}}^{N_{\mathrm{dep}}}\mathrm{\Phi }\left({\displaystyle \frac{H_k+\mathrm{\Delta }H-d_{\mathrm{1}}\cdot E_{\mathrm{sim},i}}{{\mathit{\sigma }}_{\mathrm{res}}\left(d_{\mathrm{1}}\cdot E_{\mathrm{sim},i}\right)}}\right)\\ & -\mathrm{\Phi }\left({\displaystyle \frac{H_k-\mathrm{\Delta }H-d_{\mathrm{1}}\cdot E_{\mathrm{sim},i}}{{\mathit{\sigma }}_{\mathrm{res}}\left(d_{\mathrm{1}}\cdot E_{\mathrm{sim},i}\right)}}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{pr}}}}\cdot \sum _{i=\mathrm{1}}^{N_{\mathrm{dep}}}{c}_{\mathrm{sim},ik}\end{array}\end{array}$$

where φ(⋅) and Φ(⋅) are the standard normal probability density and the cumulative distribution function, respectively. The generated spectrum is normalized by the number of simulated primaries N_pr, i. e. the number of simulated primary radionuclide decays, resulting in the same physical unit of counts per second [cps] per source activity [Bq] as for the experimental net count rate spectra (cf. Sect. 2.1.2).

As a last step, a lower level discriminator (LLD) model is applied to the generated count rate spectrum $c_{\mathrm{sim},k}^{*}$ to account for the lower spectral threshold in the experimental spectra induced by the MCA:

$$\begin{array}{}\text{(7)}& c_{\mathrm{sim},k}=c_{\mathrm{sim},k}^{*}\cdot \mathrm{\Phi }\left({\displaystyle \frac{H_k-{\mathit{\mu }}_{\mathrm{LLD}}}{{\mathit{\sigma }}_{\mathrm{LLD}}}}\right)\end{array}$$

where μ_LLD and σ_LLD are LLD model specific parameters. More information on the motivation and derivation of this LLD model can be found in Appendix B. For additional information on the uncertainty analysis for the simulated net count rate spectra c_sim,k incorporating statistical and systematic estimates, the reader is referred to Appendix A.

3 Results and discussion

3.1 Spectral comparison

Fig. 3 shows the simulated and measured spectral detector response for the seven calibration sources ⁵⁷Co, ⁶⁰Co, ⁸⁸Y, ¹⁰⁹Cd, ¹³³Ba and ¹³⁷Cs using the count rate spectra summed over the four crystals.

In general, the agreement in both absolute count rate and shape is striking considering the significant deviations previous studies have reported, especially in the lower part of the spectra (Allyson and Sanderson, 1998; Kulisek et al., 2018). Moreover, the spectral range covered by the evaluated calibration sources is significantly larger than in previous studies. For example, Allyson and Sanderson (1998) and Torii et al. (2013) used only one radionuclide, i. e. ¹³⁷Cs, to validate the spectral detector response in their studies.

Despite these improvements, there are some systematic deviations between the measured and simulated spectra for certain sources, specifically around the Compton edges (CE), the backscatter peaks (BSP) and the characteristic X-ray peaks (XRP) highlighted in Fig. 3. These deviations require further discussion.

Figure 3 Spectral comparison of the experimental and simulated detector response for different radiation sources using the summed scintillation crystal spectra. (a) ⁵⁷Co. (b) ⁶⁰Co. (c) ⁸⁸Y. (d) ¹⁰⁹Cd. (e) ¹³³Ba. (f) ¹³⁷Cs. (g) ¹⁵²Eu. Uncertainties are provided as 1 SD values (σ_tot,exp, σ_tot,sim) for the experimental as well as the simulated spectra (only visible in Fig. 3b). Backscatter peaks (BSP), characteristic X-ray peaks (XRP) and Compton edges (CE) are highlighted.

3.2 Statistical analysis

To quantify the relative spectral error over the entire spectra between simulations and measurements, the corresponding adjusted box plots (Hubert and Vandervieren, 2008) shown in Fig. 4a were computed for all calibration sources using the Library for Robust Analysis (LIBRA) (Verboven and Hubert, 2005). The spectral range considered for the adjusted box plots was defined as ${\mathit{\mu }}_{\mathrm{LLD}}+\mathrm{2}\cdot {\mathit{\sigma }}_{\mathrm{LLD}}\le H_k\le d_{\mathrm{1}}\cdot E_{\mathrm{FEP},\max }+\mathrm{2}\cdot {\mathit{\sigma }}_{\mathrm{res}}$. For E_FEP,max, the maximum photon energy of the corresponding source using a detection threshold of 5 % with respect to the emission probability of the strongest line within the associated decay branch was used (cf. Table A1 in Appendix A).

Figure 4 Model accuracy and precision analysis for the summed scintillation crystal spectra. (a) Adjusted box plot characterizing the statistical distribution of the relative spectral error between the experimental and simulated spectral detector response for various radiation sources. The statistical analysis was performed using the Library for Robust Analysis (LIBRA) (Verboven and Hubert, 2005). (b) Uncertainty quantification for the ⁶⁰Co spectral detector response. The mean count rates (c_exp, c_sim) together with the corresponding uncertainty estimates (σ_tot,exp, σ_stat,sim, σ_sys,sim) using 1 SD values are shown both for the experimental and simulated spectral detector response. Distinct spectral regions, i. e. the backscatter peaks (BSP), the Compton edges (CE) as well as the full energy peaks (FEP) are highlighted. The normalized residual level $|c_{\text{exp}}-c_{\text{sim}}|/{\mathit{\sigma }}_{\text{tot}}$ with ${\mathit{\sigma }}_{\text{tot}}:=\sqrt{{\mathit{\sigma }}_{\text{tot,exp}}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{tot,sim}}^{\mathrm{2}}}$ for a coverage factor of 2 is marked with the horizontal black dash-dotted line.

This statistical analysis shows a median relative spectral error < 10 % for five of the seven sources with respect to the defined spectral range. In contrast, the two low energy photon emitting sources ⁵⁷Co and ¹⁰⁹Cd exhibit a considerably larger median relative spectral error of about 12 % and 27 %, respectively.

This increased deviation can be explained by two factors. First, the empirical detector models, i. e. the energy calibration in Eq. (3) and detector resolution in Eq. (4), are based on calibration sources with the smallest evaluated FEP at 88.0336(1) keV. Consequently, the accuracy of the derived models deteriorates at lower spectral energies and the systematic uncertainty increases, especially due to discontinuities in the calibration and resolution models as a result of the characteristic X-ray lines of the NaI(Tl) itself (Bissaldi et al., 2009). This might also be the reason for the systematic deviations in the XRPs at the very low end of the count rate spectra highlighted in Fig. 3. Second, for the ¹⁰⁹Cd source, a custom user routine had to be used (cf. Sect. 2.2.1), which only incorporates photon but no electron emission and related bremsstrahlung. This simplification deteriorates the accuracy and therefore increases the relative spectral error even further.

We want to point out that some of the previous studies adopted integrated detector channels around the FEP to estimate the model accuracy (Allyson and Sanderson, 1998; Billings and Hovgaard, 1999). This error metric is expected to give smaller uncertainty estimates compared to channel-resolved error metrics because the goodness-of-fit is only evaluated around the FEP neglecting most of the remaining spectrum dominated by scattered photons. For example, using the FEP integrated error metric, the relative spectral error for ⁸⁸Y is equal to 1 % compared to the median relative spectral error of 5 % shown in the box plot in Fig. 4a.

3.3 Uncertainty analysis

To quantitatively assess the statistical significance of the deviations in the BSP and the CE, a detailed spectral uncertainty analysis is presented in Fig. 4b for the representative ⁶⁰Co source, which includes both features of interest not obscured by other spectral features. Multiple interesting inferences can be drawn from this analysis.

First, the systematic uncertainty for the simulated spectrum c_sys,sim is the dominant contributor to the total uncertainty in the FEP and the BSP exceeding the statistical uncertainty c_stat,sim by more than 1 order of magnitude. This confirms the conjecture of this study that a thorough uncertainty analysis including systematic uncertainties is essential for Monte Carlo based gamma-ray response simulations of AGRS systems.

Second, the residuals normalized by the total uncertainty ${\mathit{\sigma }}_{\mathrm{tot}}:=\sqrt{{\mathit{\sigma }}_{\mathrm{tot},\exp }^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{tot},\mathrm{sim}}^{\mathrm{2}}}$ and presented below the main spectral plot in Fig. 4b are consistently smaller than 2 indicated by the horizontal dash-dotted line except for the BSP, the CE and the very low part of the spectrum below about 100 keV. This implies that the difference between the simulated and measured FEP as well as the major part of the Compton continuum are well within the statistical and systematic uncertainties of the developed model and measurements considering a coverage factor of 2.

In contrast, the difference between the simulated and measured CE and BSP regions are statistically significant. It is worth adding that these findings are also supported by the statistical analysis of the ⁸⁸Y and ¹³⁷Cs radiation sources. The reader is referred to Appendix A for more details on these results as well as on the computation of the statistical and systematic uncertainties for both experimental and simulated count rate spectra.

3.4 Sensitivity analysis

As discussed in Sect. 2.2.1, besides the high-fidelity reference model (model A) used for all main simulations, five simplified models were defined to assess the sensitivity of specific model elements on the spectral detector response (model B–F). These additional models were evaluated for the representative ⁶⁰Co source and the resulting spectra are presented together with the reference model and the measured spectrum in Fig. 5d. To quantify the relative spectral error between these different models and the measured spectrum, adjusted box plots were computed for three different spectral windows highlighting the low, middle and high part of the corresponding spectrum. These box plots are shown in Fig. 5a–c together with the spectral windows highlighted in Fig. 5d.

Figure 5 Monte Carlo model sensitivity analysis for the summed ⁶⁰Co scintillation spectra. (a)–(c) Adjusted box plots characterizing the statistical distribution of the relative spectral error between the experimental and simulated spectral detector response for different Monte Carlo models (A: reference; B: no laboratory elements; C: no detector box elements; D: no source elements; E: no additional laboratory equipment elements; F: MgO reflector) and three spectral windows: (a) low: ${\mathit{\mu }}_{\mathrm{LLD}}+\mathrm{2}\cdot {\mathit{\sigma }}_{\mathrm{LLD}}\le H_k\le d_{\mathrm{1}}\cdot \mathrm{300}\mathrm{keV}$, middle: (b) $d_{\mathrm{1}}\cdot \mathrm{300}\mathrm{keV}\le H_k\le d_{\mathrm{1}}\cdot \mathrm{1000}\mathrm{keV}$ and (c) high: $d_{\mathrm{1}}\cdot \mathrm{1000}\mathrm{keV}\le H_k\le d_{\mathrm{1}}\cdot E_{\mathrm{FEP},\max }+\mathrm{2}\cdot {\mathit{\sigma }}_{\mathrm{res}}$. The statistical analysis was performed using the Library for Robust Analysis (LIBRA) (Verboven and Hubert, 2005). (d) Mean count rate spectra for the ⁶⁰Co radiation measurement (c_exp) and the different Monte Carlo models ($c_{\text{sim}}^{\text{A–F}}$) described above. The three spectral windows used to evaluate the box plots in Fig. 5a–b are highlighted.

The results from these sensitivity analysis can be summarized as follows: In the low spectral window featuring the backscatter peak and tail of the Compton continuum, the laboratory room (model B) has by far the highest impact on the spectrum followed by the detector box (model C) and source materials (model D). In the middle spectral window featuring the Compton continuum, the detector box (model C) shows the highest sensitivity followed by the source materials (model D). For the high spectral window containing the two FEPs, only the detector box materials (model C) seem to have a significant sensitivity on the spectral detector response. The two remaining models, i. e. the laboratory equipment (model E) and the reflector (model F) have negligible impact on the detector response for all three spectral windows.

These findings imply that the laboratory room, the source and the detector box elements are essential parts of the Monte Carlo model for an accurate detector response simulation. In contrast, the additional laboratory equipment and the reflector seem to be of less importance considering the adopted source-detector configuration. It is worth noting that the calibration sources in this study exhibit comparably low attenuation and self-shielding. Consequently, for sources with metal shields or bigger active volume sources, the sensitivity of the source elements on the spectral detector response is expected to be even higher.

3.5 Scintillator non-proportionality

Based on this sensitivity analysis results from the previous section, the bias in the BSP can be explained by the systematic uncertainties in the laboratory room, detector box and source elements. On the other hand, the CE region seems to be unaffected by all analyzed model elements (cf. Fig. 5d).

Interestingly, related studies adopting specialized Monte Carlo codes for NaI(Tl) scintillation detector response simulations under laboratory conditions have reported exactly the same deviations in the CE (Saito and Moriuchi, 1981; Shi et al., 2002; Gardner and Sood, 2004). Using electron spectroscopy and Compton coincidence measurements, this bias in the CE region could be attributed to the scintillator non-proportionality of NaI(Tl) and corrected for by implementing non-proportionality light yield models for the electron response in the scintillation crystals (Zerby et al., 1961; Saito and Moriuchi, 1981; Cano-Ott et al., 1999; Shi et al., 2002; Gardner and Sood, 2004; Payne et al., 2009; Moses et al., 2012). In addition, it was found that the deterioration of pulse-height spectra by the scintillator non-proportionality increases with the volume of the scintillation crystals due to the increased number of ways that the photons can interact with the crystal (Cano-Ott et al., 1999).

These findings suggest that the neglected scintillator non-proportionality is also in this study the cause for the observed bias in the CE. This conjecture is supported by the fact that standard AGRS systems contain typically the largest commercially available prismatic scintillation crystals and based on the results of other studies (Cano-Ott et al., 1999), a significant effect of the scintillator non-proportionality on the spectral detector response is expected for these systems.

4 Conclusions

The Monte Carlo model developed to simulate the spectral detector response of a standard AGRS system shows excellent accuracy and precision with a median relative spectral error < 10 % for most of the radionuclides, i. e. ⁶⁰Co, ⁸⁸Y, ¹³³Ba, ¹³⁷Cs and ¹⁵²Eu. This is a significant improvement over previously developed AGRS detector models, which possessed a relative spectral error > 10 % (Allyson and Sanderson, 1998). This superiority in model accuracy and precision can be explained by the high-fidelity physics models adopted in the FLUKA code (Battistoni et al., 2015) compared to the simplified physics models applied in the in-house codes used in previous studies (Allyson and Sanderson, 1998; Billings and Hovgaard, 1999) as well as improved accuracy in the geometric and material properties of the corresponding Monte Carlo model. This implies that, due to the significant improvement in the radiation transport code capabilities, e. g. for FLUKA, EGSnrc, GEANT or MCNP, over the past two decades, the in-house codes specifically developed for AGRS systems (Allyson and Sanderson, 1998; Billings and Hovgaard, 1999) became to a certain degree obsolete.

Moreover, compared to previous studies (Allyson and Sanderson, 1998; Billings and Hovgaard, 1999; Kulisek et al., 2018), a significant improvement in the model accuracy for the low spectral range  < 100 keV was achieved by incorporating the surrounding laboratory elements together with a detector specific LLD model. Consequently, the developed model can also be adopted to calibrate AGRS systems for low photon energy emitting radionuclides such as ²⁴¹Am or ¹⁰⁹Cd, which was not possible with previously developed models (Kulisek et al., 2018).

A sensitivity analysis confirmed the conjecture that the laboratory room, the source and the detector box elements are essential parts of the Monte Carlo model for an accurate detector response simulation. In consequence, neglecting the detector box elements, i. e. equipment surrounding the scintillation crystals such as aluminum casings, PMT, MCA or insulating foam as it was done by Allyson and Sanderson (1998) and Billings and Hovgaard (1999), seems to be questionable at best.

Based on a thorough uncertainty analysis incorporating statistical and systematic uncertainties, two important conclusions could be drawn. First, the systematic uncertainty for the simulated spectrum is the dominant contributor to the total uncertainty in the FEP and the BSP exceeding the statistical uncertainty by more than 1 order of magnitude. This systematic uncertainty is likely to be even larger for previous studies due to the reduced number of calibration sources used to derive the empirical detector models in those studies. Consequently, the incorporation of systematic uncertainties is essential in the assessment of the statistical significance for detector response simulations of AGRS systems. Second, deviations between the simulations and the measurements in the FEP as well as in the major part of the Compton continuum are statistically insignificant.

In contrast, differences in the BSP and the CE regions are statistically significant and can be attributed to the systematic uncertainties in the laboratory room model on the one hand and the neglected scintillator non-proportionality on the other hand. Previous studies have shown that the effect of the scintillator non-proportionality on the spectral detector response is increasing for larger scintillation crystals (Cano-Ott et al., 1999). These findings imply that the linear energy deposition model, which was the standard for previously developed simulation models (Allyson and Sanderson, 1998; Billings and Hovgaard, 1999; Torii et al., 2013; Sinclair et al., 2016; Kulisek et al., 2018; Zhang et al., 2018), is not appropriate for standard AGRS systems with 10.2 cm × 10.2 cm × 40.6 cm prismatic NaI(Tl) scintillation crystals in general.

The developed and validated detector model presented herein is an important step towards a simulation-based calibration methodology for AGRS systems. In contrast to previous studies (Billings and Hovgaard, 1999; Torii et al., 2013; Sinclair et al., 2016; Zhang et al., 2018), our detector model allows not only FEP characterization but also accurate simulation-based spectral predictions over the entire spectral domain relevant for AGRS applications. This prediction power can be utilized to build a detector response matrix as a function of the spectral energy and relative angular direction of the gamma-ray photons. In combination with an aircraft and an environmental model, standard spectra could be computed and used for full spectrum analysis (FSA) of the measured spectral data (Dickson et al., 1981; Grasty et al., 1985; Minty et al., 1998). Compared to the standard window evaluation methodology (Erdi-Krausz et al., 2003), FSA is particularly promising for low count rate applications like AGRS as not only the events in the FEPs are used but all the available spectral information is processed. Consequently, improved accuracy and precision for the identification and quantification of terrestrial radionuclide sources as well as reduced minimum detectable activity values can be expected. Moreover, thanks to the simulation-based approach, the standard AGRS systems could be calibrated for a significantly extended range of radionuclides such as ²⁴¹Am, ¹³¹I or ¹⁹²Ir. This extension is of high relevance for homeland security applications against nuclear terrorism as well as radiation protection applications in case of orphan source incidents or severe nuclear accidents (Acton et al., 2007; Torii et al., 2013).

In summary, the developed model shows superior accuracy and precision compared to previous models and is an important step towards a simulation-based calibration of AGRS systems. Moreover, thorough model validation and uncertainty analysis revealed statistically significant deviations in the BSP and the CE regions. Embedding of scintillator non-proportionality models in the Monte Carlo codes, characterizing and validating the angular spectral gamma-ray response as well as incorporating material and geometrical uncertainties in the systematic uncertainty analysis are necessary steps to improve the model accuracy and uncertainty estimates further and correct for the discussed biases.

Appendix A: Uncertainty analysis

Following the standard methodology of uncertainty and uncertainty propagation (Abernethy et al., 1985; Knoll, 2010), the statistical and systematic uncertainties for the experimental and simulated count rate spectra were determined. For the radiation measurements, the statistical uncertainty of the net count rate spectra c_exp,k (Eq. 1) characterized by the standard deviation was computed adopting a probabilistic Poisson model (Knoll, 2010):

$$\begin{array}{}\text{(A1)}& {\mathit{\sigma }}_{\mathrm{pois},\exp ,k}=\sqrt{{\displaystyle \frac{C_{\mathrm{gr},k}}{t_{\mathrm{gr}}^{\mathrm{2}}}}+{\displaystyle \frac{C_{\mathrm{bg},k}}{t_{\mathrm{bg}}^{\mathrm{2}}}}}\end{array}$$

where C_gr,k and C_bg,k are the gross and background count spectra in channel k together with the gross and background measurement live time t_gr and t_bg, respectively. The small statistical uncertainty in the live time measurement is neglected. The uncertainty induced by the source activity A normalization in Eq. (1) is quantified using the standard error propagation methodology for independent variables (Abernethy et al., 1985):

$$\begin{array}{}\text{(A2)}& {\mathit{\sigma }}_A=A\cdot \sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{A_{\mathrm{0}}}}{A_{\mathrm{0}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\log \left(\mathrm{2}\right)\cdot \mathrm{\Delta }t}{t_{\mathrm{1}/\mathrm{2}}^{\mathrm{2}}}}\cdot {\mathit{\sigma }}_{t_{\mathrm{1}/\mathrm{2}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\log \left(\mathrm{2}\right)}{t_{\mathrm{1}/\mathrm{2}}}}\cdot {\mathit{\sigma }}_{\mathrm{\Delta }t}\right)}^{\mathrm{2}}}\end{array}$$

with the reference activity A₀ and associated uncertainty ${\mathit{\sigma }}_{A_{\mathrm{0}}}$ provided by the vendor, the source half life $t_{\mathrm{1}/\mathrm{2}}$ and associated uncertainty ${\mathit{\sigma }}_{t_{\mathrm{1}/\mathrm{2}}}$ provided by Pearce (2008) (cf. Table A1) as well as the time difference $\mathrm{\Delta }t=t-t_{\mathrm{0}}$ between the reference date t₀ and the measurement date t together with the associated uncertainty σ_Δt. To compute the source activity A as a function of the measurement date t, the fundamental exponential law of decay is adopted, i. e. $A=A_{\mathrm{0}}\cdot {\mathrm{2}}^{-\frac{\mathrm{\Delta }t}{t_{\mathrm{1}/\mathrm{2}}}}$ (Knoll, 2010). The total experimental uncertainty is summarized as follows according to Abernethy et al. (1985):

$$\begin{array}{}\text{(A3)}& {\mathit{\sigma }}_{\mathrm{tot},\exp ,k}=\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{pois},\exp ,k}}{A}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{c_{\exp ,k}}{A}}\cdot {\mathit{\sigma }}_A\right)}^{\mathrm{2}}}\end{array}$$

For the simulations, the statistical uncertainty of the net count rate spectrum c_sim,k (cf. Eq. 7) was computed adopting the sample standard deviation (Knoll, 2010) as follows:

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& {\mathit{\sigma }}_{\mathrm{stat},\mathrm{sim},k}=\\ & \sqrt{{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{pr}}\left(N_{\mathrm{pr}}-\mathrm{1}\right)}}\cdot \left[\left(N_{\mathrm{pr}}-N_{\mathrm{dep}}\right)\cdot c_{\mathrm{sim},k}^{\mathrm{2}}+\sum _{i=\mathrm{1}}^{N_{\mathrm{dep}}}{\left(c_{\mathrm{sim},ik}-c_{\mathrm{sim},k}\right)}^{\mathrm{2}}\right]}\end{array}\end{array}$$

where c_sim,ik are the individual broadened energy deposition events in the detector channel k (cf. Eq. 6), N_dep the number of recorded events in a specific crystal and N_pr the number of simulated primaries. It is good practice in Monte Carlo studies to report not only the estimated uncertainty in the sample mean c_sim,k using the sample standard deviation ${\mathit{\sigma }}_{\mathrm{stat},\mathrm{sim},k}$ but also the so called variance of the sample variance VOV to quantify the statistical uncertainty in ${\mathit{\sigma }}_{\mathrm{stat},\mathrm{sim},k}$ itself (Forster et al., 1994):

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}{\mathrm{VOV}}_k& ={\displaystyle \frac{\mathrm{Var}\left({\mathit{\sigma }}_{\mathrm{stat},\mathrm{sim},k}^{\mathrm{2}}\right)}{{\mathit{\sigma }}_{\mathrm{stat},\mathrm{sim},k}^{\mathrm{4}}}}\\ & ={\displaystyle \frac{\left(N_{\mathrm{pr}}-N_{\mathrm{dep}}\right)\cdot c_{\mathrm{sim},k}^{\mathrm{4}}+{\sum }_{i=\mathrm{1}}^{N_{\mathrm{dep}}}{\left(c_{\mathrm{sim},ik}-c_{\mathrm{sim},k}\right)}^{\mathrm{4}}}{{\left[\left(N_{\mathrm{pr}}-N_{\mathrm{dep}}\right)\cdot c_{\mathrm{sim},k}^{\mathrm{2}}+{\sum }_{i=\mathrm{1}}^{N_{\mathrm{dep}}}{\left(c_{\mathrm{sim},ik}-c_{\mathrm{sim},k}\right)}^{\mathrm{2}}\right]}^{\mathrm{2}}}}\\ & -{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{pr}}}}\end{array}\end{array}$$

Table A1 Physical properties of the calibrated radiation sources adopted for detector characterization and simulation model validation measurements. Uncertainty values are stated as 1 SD values and provided as last significant figures, e. g. 0.12(3) corresponds to 0.12 ± 0.03.

^a Half life values are from Pearce (2008). ^b The activity of the calibrated radiation sources is computed for the corresponding measurement date based on the calibration certificate. ^c The photon energies are provided for a minimum energy of 50 keV and a detection threshold of 5 % with respect to the emission probability of the strongest line within the corresponding decay branch using the ENDF/B-VIII.0 nuclear data file library (Brown et al., 2018).

The propagation of the systematic uncertainties in the post-processing steps discussed in Sect. 2.2.2 was determined with a Monte Carlo sampling technique. The probabilistic model for the input parameters, i. e. the empirical model parameters interpreted as a random vector $\mathit{X}=(D_{\mathrm{1}},B_{\mathrm{1}},B_{\mathrm{2}}{)}^{\top }\in {\mathbb{R}}^{\mathrm{3}}$ (cf. Eqs. 3 and 4), is defined as a multivariate normal distribution $\mathit{X}\sim \mathcal{N}(\mathit{\mu }=(d_{\mathrm{1}},b_{\mathrm{1}},b_{\mathrm{2}}{)}^{\top },{\mathit{\sigma }}^{\mathrm{2}}=({\mathit{\sigma }}_{d_{\mathrm{1}}}^{\mathrm{2}},{\mathit{\sigma }}_{b_{\mathrm{1}}}^{\mathrm{2}},{\mathit{\sigma }}_{b_{\mathrm{2}}}^{\mathrm{2}}{)}^{\top })$ adopting the mean and standard deviations derived in the regression analysis (cf. Sect. 2.1.3). The $N_{\mathrm{MC}}\in {\mathbb{N}}_{>\mathrm{1}}$ independently drawn input samples $\mathcal{X}=({\mathit{x}}^{\left(\mathrm{1}\right)},\mathrm{\dots },{\mathit{x}}^{\left(m\right)},\mathrm{\dots }{\mathit{x}}^{\left(N_{\mathrm{MC}}\right)}{)}^{\top }$ are then propagated through the postprocessing pipeline, i. e. the three steps described in Sect. 2.2.2, to obtain the corresponding spectral count rate samples $\mathcal{Y}=(c_{\mathrm{sim},k}^{\left(\mathrm{1}\right)},\mathrm{\dots },c_{\mathrm{sim},k}^{\left(m\right)},\mathrm{\dots },c_{\mathrm{sim},k}^{\left(N_{\mathrm{MC}}\right)}{)}^{\top }$ with $k\in \mathit{\{}\mathrm{1},\mathrm{\dots },\mathrm{1024}\mathit{\}}$. These samples can then be used to compute the sample standard deviation ${\mathit{\sigma }}_{\mathrm{sys},\mathrm{sim},k}$ and by doing so quantify the systematic uncertainty with respect to the empirical model parameters d₁, b₁ and b₂. The total uncertainty characterized by the sample standard deviation can be summarized in the same way as for the experimental uncertainty according to Abernethy et al. (1985):

$$\begin{array}{}\text{(A6)}& {\mathit{\sigma }}_{\mathrm{tot},\mathrm{sim},k}=\sqrt{{\mathit{\sigma }}_{\mathrm{stat},\mathrm{sim},k}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{sys},\mathrm{sim},k}^{\mathrm{2}}}\end{array}$$

We want to point out that, given the spectral trends in the sensitivity analysis presented in Sect. 3.4, the qualitative conclusions are not expected to be overshadowed by the neglected systematic uncertainties in this study.

In Fig. A1, the uncertainty analysis discussed in Sect. 3.3 was repeated for other sources, which, similar to ⁶⁰Co, contain both spectral features of interest, i. e. the backscatter peak (BSP) as well as the Compton edge (CE), not obscured by other spectral features.

Figure A1 Uncertainty quantification of the spectral detector response for the summed scintillation crystal spectra. (a) ¹³⁷Cs radiation source. (b) ⁸⁸Y radiation source. The mean count rate spectra (c_exp, c_sim) together with the corresponding uncertainty estimates (σ_tot,exp, σ_stat,sim, σ_sys,sim) using 1 SD values are shown both for the experimental and simulated spectral detector response. Distinct spectral regions, i. e. the backscatter peaks (BSP), the Compton edges (CE) as well as the full energy peaks (FEP) are highlighted. The normalized residual level $|c_{\text{exp}}-c_{\text{sim}}|/{\mathit{\sigma }}_{\text{tot}}$ with ${\mathit{\sigma }}_{\text{tot}}:=\sqrt{{\mathit{\sigma }}_{\text{tot,exp}}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{tot,sim}}^{\mathrm{2}}}$ for a coverage factor of 2 is marked with the horizontal black dash-dotted line.

Appendix B: Lower level discriminator model

As already reported in previous studies (Kulisek et al., 2018), at the very low end of the spectral detector range, typically below 100 keV for AGRS systems, the experimental spectra are significantly lower than the model predictions. This discrepancy can be mainly attributed to the lower level discriminator (LLD) adopted in the MCA (Knoll, 2010). So far, no attempts have been found in the literature to include this effect in a consistent way for AGRS Monte Carlo models. In this study, we apply a Gaussian model to characterize this lower level discriminator:

$$\begin{array}{}\text{(B1)}& c_{\mathrm{sim},k}=c_{\mathrm{sim},k}^{*}\cdot \mathrm{\Phi }\left({\displaystyle \frac{H_k-{\mathit{\mu }}_{\mathrm{LLD}}}{{\mathit{\sigma }}_{\mathrm{LLD}}}}\right)\end{array}$$

where Φ(⋅) is the standard normal cumulative distribution function and μ_LLD, σ_LLD are the LLD model parameters. This heuristic model is motivated by the idea that the fixed lower level threshold is fluctuating due to inherent electronic noise in the MCA and associated electronics resulting in a normally distributed cut-off with the mean equal to the predefined LLD threshold and a standard deviation determined by the level of electronic noise.

Figure B1 ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}$ contour maps for the four individual scintillation crystals 1–4 (a–d) as a function of the lower level discriminator model parameters μ_LLD and σ_LLD . The global minima (${\mathit{\chi }}_{\mathit{\nu },\min }^{\mathrm{2}}$) are marked with a plus sign together with the corresponding 1 and 2-σ confidence regions indicated by the dashed-dotted and dashed contour lines, respectively.

The LLD model parameters μ_LLD, σ_LLD are determined by a numerical optimization approach, i. e. simulated spectra, which include the LLD model, are compared with experimental spectra and the model parameters are optimized according to some objective function. In this study, we adopt the Chi-squared statistic (${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}$) as the objective function to solve the global minimization problem:

$$\begin{array}{}\text{(B2)}& \left({\widehat{\mathit{\mu }}}_{\mathrm{LLD}},{\widehat{\mathit{\sigma }}}_{\mathrm{LLD}}\right)=\underset{{\mathit{\mu }}_{\mathrm{LLD}},{\mathit{\sigma }}_{\mathrm{LLD}}}{\arg \min }{\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{LLD}},{\mathit{\sigma }}_{\mathrm{LLD}}\right)\end{array}$$

where the ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}$ objective function is defined as follows:

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& {\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\left({\mathit{\mu }}_{\mathrm{LLD}},{\mathit{\sigma }}_{\mathrm{LLD}}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{N_{k^{\prime }}-\mathrm{2}}}\sum _{k^{\prime }}={\mathrm{1}}^{N_{k^{\prime }}}{\displaystyle \frac{{\left(c_{\exp ,k^{\prime }}-c_{\mathrm{sim},k^{\prime }}\left({\mathit{\mu }}_{\mathrm{LLD}},{\mathit{\sigma }}_{\mathrm{LLD}}\right)\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{tot},\exp ,k^{\prime }}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{tot},\mathrm{sim},k^{\prime }}^{\mathrm{2}}}}\end{array}\end{array}$$

with $\mathit{\{}{k}^{\prime }:H_{k^{\prime }}/d_{\mathrm{1}}\le \mathrm{50}\mathrm{keV}\mathit{\}}$. The global minimization problem in Eq. (B2) is solved by a radial basis surrogate solver (Regis and Shoemaker, 2007) provided by the MATLAB^® Global Optimization Toolbox. The sample size of the solver was set to 10³ and the convergence was checked by convergence plots. The ⁶⁰Co measurement and simulation model was selected for the LLD optimization procedure because this source was the only one measured in this study that does not show any characteristic X-ray emission above 10 keV, which could obscure the LLD threshold. The resulting ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}$ contour maps for the four crystals are presented in Fig. B1. For each crystal, the global minimum (${\mathit{\chi }}_{\mathit{\nu },\min }^{\mathrm{2}}$) marked by the plus sign was successfully found. Confidence limits are illustrated adopting 1 and 2 SD contour lines (Avni, 1976). As it emerges from this analysis, the mean LLD parameter for the detector crystal 1 is significantly smaller compared to the crystals 2, 3 and 4. The spectra presented in Fig. 3 confirm that the LLD model corrects not only the LLD effect for the ⁶⁰Co source adequately but also for the other sources, which were not used to fit the LLD model parameters. For the reader's convenience, the optimized LLD model parameters ${\widehat{\mathit{\mu }}}_{\mathrm{LLD}}$ and ${\widehat{\mathit{\sigma }}}_{\mathrm{LLD}}$ will be referred to as μ_LLD and σ_LLD in the remainder of this study.

Appendix C: Nomenclature

Symbols
A	Source activity [Bq]
a	Normalizing Gaussian model parameter [cps Bq−1]
B1, B2	Random spectral resolution model variables [−]
b1, b2	Realizations of B1, B2 [−]
C	Count spectrum [−]
c	Count rate spectrum normalized by the source activity A [cps Bq−1]
D1	Random slope variable for the spectral energy calibration model [keV−1]
d1	Realization of D1 [keV−1]
E	Spectral energy [keV]
Edep	Energy deposition [keV]
Esim	Simulated energy deposition [keV]
H	Channel number [−]
k	Coverage factor [−]
μ	Mean Gaussian model parameter [−]
N	Cardinality of a set [−]
R2	Coefficient of determination [−]
σ	Standard deviation Gaussian model parameter [−]
t	Measurement live time [s]
$t_{\mathrm{1}/\mathrm{2}}$	Source half life [s]
X	Random model parameter vector
𝒳	Random model parameter set
x	Realization of X
𝒴	Random response variable set
Mathematical operators
FWHM (⋅)	Full width at half maximum
$\mathcal{N}\left(\cdot ,\cdot \right)$	Normal distribution
VOV (⋅)	Variance of the sample variance
μ(⋅)	Sample mean
σ(⋅)	Sample standard deviation
σres(⋅)	Standard deviation resolution model
Φ(⋅)	Standard normal cumulative distribution function
φ(⋅)	Standard normal probability density
${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\left(\cdot \right)$	Reduced Chi–squared statistic

Subscripts
0	Reference value
b	Experimental baseline component
bg	Experimental background component
exp	Experimental net component
gr	Experimental gross component
i	Energy deposition event index
k	Channel index
l	Full energy peak index
MC	Monte Carlo
max	Maximum
min	Minimum
pk	Gaussian peak
pois	Poisson model component
pr	Simulated primary
res	Resolution
sim	Simulated net component
stat	Statistical uncertainty component
sys	Systematic uncertainty component
tot	Total uncertainty component combining statistical and systematic uncertainty components

Superscripts
′	Subset
$\widehat{}$	Estimator
*	Component without LLD correction
m	Monte Carlo sample index
Acronyms
AGRS	Airborne gamma-ray spectrometry
BSP	Backscatter peak
CE	Compton edge
cps	Counts per second
FEP	Full energy peak
FSA	Full spectrum analysis
IAEA	International Atomic Energy Agency
LIBRA	Library for Robust Analysis
LLD	Lower level discriminator
MCA	Multi-channel analyzer
OAT	One at a time
PE	Polyethylene foam
PI	Prediction interval
PMT	Photomultiplier tube
PTFE	Polytetrafluoroethylene
SD	Standard deviation
WNLLS	Weighted non-linear least-squares
XRP	Characteristic X-ray peak