Footprint and height corrections for UAV-borne gamma-ray spectrometry studies

Advancements in the development of gamma-ray spectrometers (GRS) have led to small and lightweight spectrometers that can be used under unmanned aerial vehicles (UAVs). Airborne GRS measurements are used to determine radionuclide concentrations in the ground, among which the natural occurring radionuclides 40K, 238U, and 232Th. For successful applications of these GRS sensors, it is important that absolute values of concentrations can be measured. To extract these absolute radionuclide concentrations, airborne gamma-ray data has to be corrected for measurement height. However, the current analysis models are only valid for the height range of 50–250 m. The purpose of this study is to develop a procedure that correctly predicts the true radionuclide concentration in the ground when measuring in the UAV operating range of 0–40 m. An analytical model is developed to predict the radiation footprint as a function of height. This model is used as a tool to properly determine a source-detector geometry to be used in Monte-Carlo simulations of detector response at various elevations between 0 and 40 m. The analytical model predicts that the smallest achievable footprint at 10 m height lies between 22 and 91 m and between 40 and 140 m at 20 m height. By using Monte-Carlo simulations it is shown that the analytical model correctly predicts the reduction in full energy peak gamma-rays, but does not predict the Compton continuum of a spectrum as a function of height. Therefore, Monte-Carlo simulations should be used to predict the shape and intensity of gamma-ray spectra as a function of height. A finite set of MonteCarlo simulations at intervals of 5 m were used for the analysis of GRS measurements at heights up to 35 m. The resulting radionuclide concentrations at every height agree with the radionuclide concentration measured on


Introduction
Technological advancements in the development of unmanned aerial vehicles (UAVs) have led to the availability of affordable drones that can be used for geophysical gamma-ray spectrometry surveys. These UAVs can carry payloads in the order of kilograms, which is roughly the weight of a gamma-ray spectrometer needed to efficiently map radionuclide activities in an area (Nicolet and Erdi-Krausz, 2003). UAV based surveys combine the advantageous properties of land-based and airborne studies.
Ground-based surveys are typically used to map an area with a high spatial resolution, but are limited by the accessibility of the terrain.
Airborne surveys overcome this limitation and can conveniently be used to cover rocky, wet, and densely vegetated terrain. However, manned airborne surveys have the inherent drawback that they are expensive and have to follow the aircraft safety rules, which impose a minimum flying height and therefore, airborne surveys cannot obtain the high standard of spatial resolution which the land-borne surveys can.
UAVs, on the other hand, can fly over areas at a relatively low height, thus collecting high-resolution spatial radionuclide information while taking advantage of the possibility to fly over to access areas. Furthermore, with the use of UAVs, as they are unmanned, it is possible to measure in areas that would pose a threat to humans because of (radioactive) pollution. Most of the recent studies found in literature that describe UAV borne gamma-ray measurements use this platform for the mapping of man-made radioactive material in the environment, such as 137 Cs (Martin et al., 2016;Sanada et al., 2016;Tang et al., 2016). These studies use a relatively small detector (max 210 ml) and identify relatively high count rates (due to radioactive pollution) which are manifested as surface or point sources. The mapping of natural occurring radionuclides is significantly different from mapping radioactive material as a result from nuclear accidents because natural sources have much lower count rates and are manifested as volume sources.
There are two commonly used methods to calculate the concentrations of naturally occurring radionuclides from the recorded gamma-ray spectra: the windows method (Nicolet and Erdi-Krausz, 2003) and full spectrum analysis (Hendriks et al., 2001). The windows method only uses the counts collected in windows set around three prominent full energy peaks of the natural occurring radionuclides 40 K, the 238 U-series, and the 232 Th-series. Full spectrum analysis uses almost all collected counts in the whole spectrum for these radionuclides. Both methods rely on knowledge about the shape of the spectra. This can be obtained with calibration measurements, as is common for the windows method, or by using Monte-Carlo simulations, as is common for the full spectrum analysis approach. The construction of a Monte-Carlo model to create standard spectra for ground-based measurement systems has become a standard procedure (Van der Graaf et al., 2011).
To translate gamma-ray spectra acquired at a certain altitude to absolute radionuclide concentrations in the soil, several analytical steps have to be taken. One of these steps is the correction for height that takes into account both the attenuation of the signal, due to the layer of air between the ground and the detector, and the change in the field of view (footprint) of the detector. Detailed instructions on how to process gamma-ray spectra by using the windows method, including approximations for flight height corrections, are available (Nicolet and Erdi-Krausz, 2003).
However, these height corrections (Nicolet and Erdi-Krausz, 2003) cannot be directly used in the UAV operating range (0-40 m) due to an approximation that is only valid in the conventional airborne operating range (50-250 m). Height corrections valid for higher altitudes have been derived in earlier studies (Duval et al., 1971;Grasty et al., 1979). These studies use a modelling approach to estimate the reduction in gamma-ray intensity and aim to locate the origin of radiation for airborne gamma-ray measurements in geophysical studies. These studies have not been adopted as the standard in the processing guidelines by the International Atomic Energy Agency (IAEA), most likely because the operating range of 0-50 m was not feasible for airborne spectral surveys at the time. Despite the IAEA height corrections only being valid for above 50 m, multiple surveys flying at lower heights have been reported using these corrections (Billings et al., 2003;Kock and Samuelsson, 2011;Pfitzner et al., 2003;Šá lek et al., 2018).
The present study focusses on height corrections for gamma-ray spectrometry surveys in the UAV operating range of 0-40 m altitude. The aim is to describe the change of the gamma-ray spectrum in the 0.3-3 MeV energy range, both in intensity and in shape. Based on this description, corrections for flight height will be suggested.
This article contains three approaches to these height corrections: analytical, computational and experimental. Firstly, the analytical approach is an extension of earlier work by Duval et al. (1971) and Grasty et al. (1979) and is a description of the change in intensity of photons that directly result from nuclear decay as a function of height. Secondly, a computational approach is used in order to describe photons that have interacted either in the ground or in the air by Compton scattering and/or pair production. In this (computational) approach the analytical results form the basis for the construction of a Monte-Carlo model that predicts the spectral shape measured at a certain height.
Finally, the results of both the analytical and the computational approach are compared with an experimental study. In this study, measurements were taken with a 2000 ml CsI scintillation detector at heights 13, 17, 21, 26, 31, and 35 m altitude above an agricultural field.

Analytical ground volume estimation
In airborne gamma-ray spectrometry studies there are two competing processes that influence the amount of radiation that is detected by the gamma-ray spectrometer, namely: 1) the change of signal due to the increase of the footprint (field of view) that contributes to the gamma-ray signal when moving to more elevated heights; and 2) the decrease in signal due to attenuation by the layer of air between the ground (the source) and the detector.
The following equation (based on Lambert's Law (Ingle and Crouch, 1988)) represents the contribution to the signal in the detector originating from a volume element in the soil that emits mono-energetic gamma rays which are attenuated by the ground and air between the volume element and the detector (Duval et al., 1971): where I total is the number of detected gamma rays per second (s − 1 ), A is the surface area of the detector (m 2 ), ε is the efficiency of detector, γ is the number of mono-energetic gamma rays emitted per unit volume of source material per second (s − 1 m − 3 ), R represents the total distance between the detector and the origin of the gamma-ray (m), r a and r g are the distances the gamma-rays travel (m) in air and the ground respectively, θ is the angle of the detector with respect to the normal of the surface pointing downwards (radians), ρ a and ρ g are the densities of the air and the ground (kg m − 3 ), μ a and μ g are the mass attenuation coefficients (m 2 kg − 1 ) in the air and ground at the energy of the emitted gamma rays, respectively. In all calculations in this paper the energy dependent mass attenuation coefficients are taken from the standard reference database of the National Institute of Standards and Technology (NIST, Gaithersburg, MD, USA; Hubbell and Seltzer, 2004). The equation for the total detected intensity I total (s − 1 ) by a detector at a height h (m) above a semi-infinite homogeneous volume is found by the integration of equation (1) over the total volume of the ground. Using the substitutions r a = h/cos(θ) and r g = R − h/cos(θ) and integrating equation (1) results in: where E 2 represents the exponential integral of order 2 (Abramowitz and Stegun, 1970). Equation (2) can be used at any altitude for height corrections in airborne gamma-ray spectrometry. Several analytical approximations where the contribution of the ground has been limited to a finite volume have been published (Duval et al., 1971;Grasty et al., 1979). This finite volume is achieved by setting limits to the integral presented in equation (1). Duval et al.'s derivation is based on assuming that the maximum distance the gamma-rays travel through the ground is constantwhich omits the effect of increasing absorption by air with increasing opening angle. Grasty et al. set a limit to the maximum angle θ to study the contribution of each concentric circle of ground. The previously published analytical models do not have the aim to precisely describe the spatial origin of the radiation and neglect the precise depth of the origin of the radiation. In the section below a new method is proposed that fully accounts for the attenuation that the gamma-rays undergo in the ground and therefore, accurately describes the depth distribution of radiation. For this new approach the integration limits of R in equation (1) are set to an isoline that can be interpreted as the line that defines the points where the total amount of attenuation in the ground and in the air is equal. This line is given by the equation: μ g ρ g r g + μ a ρ a r a = d μ g ρ g + h μ a ρ a where d is the distance in the ground directly below the detector. The dependency on θ enters in equation (3) through r g and r a in the equation can be rewritten into an equation for the integration limit of R: This line is schematically shown in Fig. 1a by the solid black isoline (R) in the ground. Using the limits ϕ : 0 →2π, R : R isoline → ∞, θ : 0→ θ 1 in the integration of equation (1) we find: Equation (5) is a new description of the intensity that uses the isoline as described by equation (3) as a boundary. Equation (5) describes the part of the volume that lies outside the isoline and has a maximum angle θ 1 . The area represented by this equation is schematically shown by the shaded area in the right part of Fig. 1a. As the distance d in the ground decreases towards zero for large angles, in our approach a maximum angle θ max is introduced. This maximum angle is defined by the point where the attenuation through the air is equal to the attenuation through the ground and air directly beneath the detector. The maximum distance through the air, r a,max follows from equation (3) by setting r g to zero: and thus, the maximum angle θ max is defined as: To calculate the intensity of the gamma signal that originates from the volume that falls within the isoline defined by equation (3) we can subtract the areas outside this isoline from the total intensity given by equation (2). The shaded area in Fig. 1a is defined by equation (5). The white area in Fig. 1a below the surface outside the isoline with θ > θ 1 has previously been described by Duval (Duval et al., 1971) and is given by: This results in the intensity that falls within the isoline and (θ < θ max ) to be described by: and the relative amount of radiation originating from this volume is given by: Equations (9) and (10) are schematically represented by the shaded area in Fig. 1b. The volume of this ground section is given by: Fig. 1. Schematic cross sections of the ground volume included in the analytical approach. The schematic representations are rotationally symmetric around the height-axis that coincides with the center of the detector. a) The shaded area is described by equation (5) for: θ max (left side) and θ 1 (right side). b) The shaded area is described by equation (9) (intensity) and equation (11) (volume). The distance R total is updated to the extent that the maximum gamma-ray attenuation remains equal. This distance in the ground r g is equal to d directly below the detector, and decreases to zero at the maximum angle θ max determined by equation (7).
where g = (d + fh) ⋅and⋅f = μ a ρ a μ g ρ g Comprehensive derivations of the equations presented in this section can be found in the supplementary material.

Ground volume estimation
The analytical method described in the previous section can be used to estimate the response of a detector when measuring above a homogenous volume. However, this analytical model can only calculate the reduction in intensity of monoenergetic peaks and does not account for photons that have interacted by Compton scattering or pair production in the soil or the air. To account for these effects, Monte-Carlo simulations were used to calculate the response of the detector to the radiation from the radionuclides 40 K, 238 U-series and 232 Th-series naturally occurring in soil. Monte-Carlo simulations are a class of computational algorithms that rely on the random sampling of events. For these simulations a finite geometry has to be defined from which the gamma-rays originate. These models have to include a sufficiently large volume of the ground to correctly model the true response. On the other hand, it is beneficial to minimize the size of the model (ground volume considered as source) to reduce calculation time. Consequently, the integrated ground volume following from the new model described in section 2.1 is estimated and used for the Monte-Carlo simulation.
For the construction of the Monte-Carlo model two parameters are of interest: the maximum depth and the radius of the radioactive source used in the model. The depth distribution in soil follows from these parameters by using equation (10).
The radius of the source follows from θ. Equations (9) and (11) are inherently limited by a maximum angle θ max that results from a depth d. This θ max results in a gradually decreasing top layer depth with increasing angle, but this top layer has a steep increase in terms of volume. For the practical implementation of the Monte-Carlo model, a cutoff angle has to be determined that bounds the geometry of the simulation. This cutoff angle has to meet the following requirements: θ cutoff ≤ θ max and the amount of radiation emitted by the volume outside θ cutoff should be insignificant.
Both the depth and maximum angle are studied as a function of height by using the relative amount of radiation that is emitted within a certain volume, as described by equation (10), versus the volume of that ground section as described by equation (11). These results will be used to determine the extent of the geometry needed for the Monte-Carlo simulations.

Monte-Carlo simulations
Full spectrum analysis (Hendriks et al., 2001) uses standard spectra to calculate the radionuclide concentrations from a measured spectrum. These standard spectra are constructed by using Monte-Carlo simulations and validated by measurements in a calibration facility (Van der Graaf et al., 2011).
Based on the well-known interactions of gamma-rays with matter, Monte-Carlo simulations can be used to estimate the transport of gamma-rays from a source (e.g. the ground) to a detector in a relevant geometry (e.g. detector at a certain height above the ground). In Monte-Carlo simulations a tradeoff has to be made between calculation time and uncertainty of the results. Simulations to calculate the response of a gamma-ray detector in ground-based geophysical situations typically take less than 24 h on a current day computer to converge to an acceptable result. Separate simulations have to be done for each of the three naturally occurring radionuclides. However, for airborne spectra, where the detector is much further away from the soil, the relevant source of the Monte-Carlo model increases significantly in volume. Computation time increases rapidly with increasing height because the soil volume that the detector sees increases roughly with the height squared. E.g. an increase from 0.80 m height to 20 m will increase the volume by a factor 625, resulting in a similar increase in computation time.
In general, when constructing a Monte-Carlo model, there is usually a tradeoff between spatial size and computational time. On the one hand the geometry of the model should be large enough so that it will correctly describe the resulting spectra. But on the other hand, the model should not be too large to avoid computation time which takes too long. An obvious approach to try and shorten computation time is by eliminating parts of the geometry that do not significantly contribute to the signal in the detector. By implementing such balanced Monte-Carlo models the change in spectral shape as a function of height can be studied with acceptable computation times.
The uncertainty in a Monte-Carlo simulation decreases with the number of radioactive decays simulated that actually lead to an event inside the detector. Since the number of decays that can be simulated in unit time is roughly independent of the model size, an increase in the size of the model either in source volume or detector-source distance will lead to a decrease in the number of photons that reach the detector in unit time. It is therefore evident that a proper limitation of the source volume is of crucial importance for effective simulations of the spatially large source-detector geometries at play in airborne surveys. Limiting the size of the source (the ground) is established as a function of detector height by requiring the relative intensity given by equation (10) to be at least 99% of the theoretical total intensity given by equation (2). The ground source was modeled as seven circular symmetric layers, where the deeper layers had smaller radii. The resulting geometry encapsulates the 99% isoline given by equation (2).
In this study MCNP 6.2 (Goorley et al., 2013) is used to determine the response of the spectrometer. MCNP is a multipurpose Monte-Carlo radiation transport code that randomly generates and tracks nearly all particles at nearly all energies. Currently, version 6.2 is the latest of a development that started at Los Alamos National Laboratory (Los Alamos, NM, USA) nearly sixty years ago. To further decrease the calculation time, various variation reduction techniques were used. The variance reduction techniques are optimizations that are embedded in the software package and the reader is directed towards the MCNP manual for an extended description of these techniques (Goorley et al., 2013).
The standard variance reduction techniques of MCNP that were implemented were: source biasing (biasing particles towards the detector), weight windows (generated by MCNP) and forced collisions in the spectrometer scintillation crystal. These variance reduction techniques apply a weight to the particles to retain proper statistics. Furthermore, an energy cutoff of 300 keV was used and the calculation was distributed over 40 computer cores.
To further decrease the simulation time the model was adapted to include a grid of detectors. Only the particles that reach a detector contribute to the result of the simulation. Therefore, increasing the number of detectors raises the number of particles that contribute to the result.
To ensure that the use of such a grid of detectors does not change the results compared to a simulation with a single detector this grid has to meet the following requirements: 1) the outer boundaries of the grid have to be small compared to the size of the model to guarantee that no geometry boundaries effects are induced; and 2) the detectors cannot change the yield of each other. This can be caused by shielding radiation, or by events caused by a detector ending up in a neighboring detector which contribute to the Compton continuum. The extent of the ground volume has appropriately been widened so that the requirement of the inclusion of 99% of the radiation is met for all detectors in the grid.
With these requirements, the results of the Monte-Carlo simulations can be divided by the number of detectors while obtaining a smaller uncertainty due to the higher number of events scored in the grid. However, implementing a grid of detectors adds additional elements to the geometry which in turn slows down the calculation speed. It was determined that the optimal number of detectors in the grid was 20 × 20 (in total 400 detectors). A further increase of this number of detectors resulted in a lower rate of scored events in the grid per unit of time.
In the simulations, the spectrometer model is based on a MS-2000 CsI spectrometer developed by Medusa Radiometrics (Medusa Radiometrics, 2019). This device contains a 2000 ml CsI scintillation crystal fitted with a photomultiplier tube, an MCA and a high voltage module all connected to the proprietary Medusa detector operating system (mDOS) where the information from the spectrometer is stored and analysed in real-time. This detector is chosen because it is used in the field measurements described in section 2.4.
The ground has been modeled as a homogenous volume of SiO 2 with a density of 1.5 g cm − 3 , mixed with 0.3 g cm − 3 H 2 O totaling to a compound density of 1.8 g cm − 3 . The density and composition of the soil is very specific for each area and typically varies within a field and even with depth. The above densities are an approximation of the complex soil composition found in a field and in line with common values reported in literature (Rabot et al., 2018).
The air has been modeled up to 40 m, or to extend at least 20 m above the detectors with a density of 1.225⋅10 − 3 g cm − 3 . Simulations have been done for heights between 5 and 40 m at an interval of 5 m. Additionally, simulations at 0.80 m height have been included because this is the typical surveying height for ground-borne gamma-ray surveys.
The requirements for using a grid of detectors were met by using a spacing of 1 m between the centers of the detectors to prevent shielding from each other. The grid was centered at a height h above the center of the ground volume. These parameters determined that the outer boundaries of the grid were given by a square of 19 m centered at the origin. The extent of the Monte-Carlo model will be determined in the next section, and one of the boundary conditions is that the spatial extent of the simulation should be large compared to this grid boundary. The ground source will be modeled such that for all detectors in the grid, at least 99% of the theoretical total intensity given by equation (2) is included (Fig. 2). To verify that the detectors do not shield each other, a simulation has been done at 0.80 m in which the response of the detectors in the outer layer of the grid are compared with the inner layer. It was found that both sets of detectors produced comparable results, and thus it can be assumed that for all heights there is no inter-detector influence.
The spectrum in the detectors was simulated by using a MCNP F8 tally that scored the energy deposition in the detector in 10 keV energy bins in the range of 0.3-3 MeV. Simulations have been done for the natural occurring radionuclide 40 K, and for the 238 U-and 232 Th-series. A simulation for 40 K was marked as sufficiently accurate if the relative error of the bin containing the 1.46 MeV peak was below 5%. 238 U-series and 232 Th-series simulations were marked sufficiently accurate if the relative errors in the bins that contained the ten most prominent peaks of the respective decay chains were below 5%.
For comparison with the field measurements the unbroadened MCNP output spectra have been broadened by using gaussian energy broad- where E represent the energy of the channel). For the parameter b a value of 0.05 is used, which is a typical for a gamma-ray spectrometer used in airborne measurements and gives a resolution of 9.4% at 662 keV. The schematic overview shows a cross section of a geometry consisting of a grid of 3 detectors. A red line in the ground (below the black horizontal line) represents the source boundaries in the ground (7 layers) used for the Monte-Carlo simulation. The left detector (#1.1 orange) has a footprint (orange line: footprint #1.1) that falls within the left source boundaries (red line), while the right detector (#1.3 blue) has a footprint (blue line: footprint #1.3) that falls within the right source boundary (red). Effectively, the lowest layer of the source is widened so that 99% of the detectors footprints fall within the source boundaries. This schematic overview only shows 3 detectors, the actual implementation uses a grid of 20 by 20 detectors. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Field measurements
Stationary measurements at several heights were taken by using a commercially available APID One UAV (Fig. 3) manufactured by MainBase (Linköping, Sweden). This platform is designed in the shape of a helicopter and uses standard petrol to power the engines. The APID One has a rotor diameter of 3.3 m, empty weight of 130 kg and a maximum take-off weight of 210 kg.
The UAV was fitted with a MS-2000 gamma-ray spectrometer. GPS and barometer were connected to the measurement system to accurately determine the position and height of the measurements.
The gamma-ray spectrometer measured continuously while the UAV hovered at 13, 17, 21, 26, 31 and, 35 m height above a fixed point in a 50 ha agricultural field on Bjertorp Farm in southwest Sweden (58.248 • N; 13.128 • E). At each height the UAV hovered for approximately 60 s. That time span was sufficient to have a statistical uncertainty of less than 5% for the radionuclide concentrations determined by using Gamman® spectral analysis software (Medusa Radiometrics, 2020) (using full spectrum analysis with 1 Bq/kg based standard spectra of 40 K, 238 U-series and 232 Th-series). The measurement footprint is height dependent and the radionuclide concentration is the average of the concentration within this footprint.
Additionally, a ground-based measurement was collected with the same MS-2000 gamma-ray spectrometer. This measurement is meant to establish the spatial homogeneity of radionuclides around the center of the UAV height measurement. This ground-based survey was conducted by mounting the spectrometer on the back of a tractor at a height of 0.80 m and measuring the radionuclide concentration up to a radius of 300 m around the center of the height measurements. The tractor mapped the field at an average speed of 5.6 m s − 1 (20 km h − 1 ), used a line spacing of 25 m and recorded spectra at 1 hz. Nuclide concentrations were established by using Gamman® with standard spectra calibrated at 0.80 m height. Spatial nuclide distribution maps were interpolated using the geostatistical method ordinary block kriging (Burrough et al., 2015).
The stationary height measurements were analysed with two different approaches. The first approach was to sum the spectra taken at the same height and analyse them using Gamman® and a Monte-Carlo generated standard spectrum at 0.80 m. The second approach was to (quadratically) interpolate Monte-Carlo generated standard spectra to generate a standard spectrum for the measurement height (standard spectra for heights 5,10,15,20,25,30, and 35 m are available).

Results
This section consists of six subsections. In subsections 3.1-3.3 the results of the analytical integration for the volume of ground that contributes to the signal in the detector are given. These results have been used to define the ground volume used in the Monte-Carlo calculations, which are presented in subsection 3.4 and compared with a scaled spectrum using equation (2) in section 3.5. The final subsection (3.6) presents the field measurements as a function of height. Throughout the results section the term relative intensity is used and will be expressed as a percentage. Fig. 4 shows the intensity with respect to the intensity at 0.80 m measured by a detector as a function of the height for the characteristic gamma energies of 40 K (1.46 MeV), the 238 U-series (1.77 MeV) and the 232 Th-series (2.62 MeV). In Fig. 5, these intensities are shown for both the IAEA model (red lines) and the model given by equation (2). The IAEA model uses an exponential function with typical experimentally determined window attenuation coefficients given in IAEA (IAEA, 1991). As expected, the reduction in intensity is stronger for lower energy gamma-rays. For the model described by equation (2), at 40 m height, the reduction is approximately 40% for the 2.62 MeV gamma of 232 Th compared to 50% for the 1.46 MeV gamma of 40 K. As described in (Nicolet and Erdi-Krausz, 2003), the IAEA model is valid for the 50-250 m range and this figure clearly shows that there is a difference in the predicted signal decrease between the IAEA and the model presented in equation (2) in the height range from 0 to 40 m.
Sections 3.4 and subsequent sections present the results of the Monte-Carlo simulations which have been done for the natural occurring radionuclides 40 K, 238 U-series and 232 Th-series. Only the results for 40 K, which emits a single gamma energy, and 232 Th, which has a decay chain are presented. It is established that the resulting reduction of the spectra that results from the 238 U-series are very similar to the reduction of the 232 Th-series, and the 232 Th-series was chosen because it contains the gamma-ray with the highest energy associated with the decay of naturally occurring radionuclides. Fig. 6 shows cross sections of the ground that contributes to the measured intensity for two detector heights. These figures use equation (10) to determine the isolines that define the cross section in the soil that correspond to relative intensities of 65%, 95% and 99%. The cross sections are rotationally symmetric and the contributing ground volume falls within the isoline. The isolines are based on a soil bulk density of Fig. 3. Photograph of the APID One UAV carrying the gamma-ray spectrometer.

Fig. 4.
Intensity with respect to 0.80 m versus height for the three characteristic energies associated with the naturally occurring radionuclides ( 40 K, 238 U, and 232 Th). The plot has been made by using equation (2). 1.8 g cm-3 (1.5 g cm-3 SiO 2 and 0.3 g cm-3 H 2 O) and the mass attenuation coefficient for the 2.62 MeV gamma (from 208 Tl in the 232 Th-series) which has the highest energy emitted by the natural occurring radionuclides. As gammas with lower energies have higher attenuation coefficients, using the 2.62 MeV gamma results in the largest ground volume contributing to the intensity measured by the detector. Fig. 6 shows that the ground volume from which the radiation gives a 99% relative intensity has a maximum depth extending not beyond 65  (2) and the IAEA model (valid for the range 50-250 m altitude) (for the three characteristic energies associated with the naturally occurring radionuclide 40 K, 238 U, and 232 Th, from left to right respectively).

Fig. 6.
Cross sections of the ground showing the origin of the radiation that gives 65%, 95% and 99% of the intensity (gamma energy 2.62 MeV, ground density of 1.8 g cm − 3 (1.5 g cm − 3 SiO 2 and 0.3 g cm − 3 H 2 O)) for two different heights above the ground. Note that the horizontal axis for figure a) extends to ±20 m and b) extends to ±450 m. The figures represent the relative contribution of the gamma-signal. The absolute signal detected as a function of height changes significantly as shown in Fig. 4. cm when the detector is placed at 0.80 m height above the surface. For larger heights, this maximum depth increases slightly to 67 cm for 20 m height and 68 cm for 40 m height for the 99% isoline. It can be seen that the isolines in the ground widen and for larger heights, the top layer of soil becomes increasingly important as a source of radiation.

Depth of Monte-Carlo simulation geometry
With regard to Fig. 6, it can be concluded that 99% of the radiation originates from within the top 70 cm of the soil for the whole altitude range of 0-40 m. Fig. 7 shows the amount of radiation as a function of the volume that is included in the analytical integration, calculated by equation (11). The three depths, 17, 45 and 70 cm have been selected because when using these maximum depths, the equation for the amount of radiation captured in the detector reaches a maximum of 65%, 95% and 99% of the intensity of a semi-infinite model, respectively. The values of the depths are constant and the volume increase is realized by increasing the angle in equation (11).
Both graphs in Fig. 7 show a steep initial increase of the relative intensity as a function of increasing volume and then flatten and converge to the maximum value for each depth. This means that the largest fraction of the detected radiation originates from the volume that is located directly below the detector. The tails of the isoline contribute only a small fraction of the total radiation intensity. Note that the volume axis has a logarithmic scale and so the tails represent a relatively large part of the volume.
From Figs. 6 and 7 it is concluded that a Monte-Carlo model that calculates the response on a ground source of 232 Th should extend down to at least 70 cm in the ground directly below the detector for the model to capture more than 99% of the maximum intensity. Fig. 6 can be used to estimate the depth and radius of the ground geometry to be used in a Monte-Carlo simulation. Fig. 6b shows that for 20 m height, the 99% isoline has a tail with a depth of still a few centimetres. Fig. 8 shows the radius of the analytical model as a function of the height. The lines in the figure represent the radius that results from solving equation (10) by setting the relative intensities to 65%, 95% and 99% and using a depth of 70 cm. The resulting angle θ can be converted to a radius by using the detector height, and the volume can be calculated by using equation (11).

Radius of Monte-Carlo simulation geometry
From Fig. 8 it is concluded that the radial coordinate of the ground volume in a Monte-Carlo model should extend to 450 m such that the relative intensity is at least 99% for heights up to 40 m. At 10 m height the 65% and 95% lines intersect at a radius of 22 and 92 m respectively. For 20 m height these points intersect at radii of 40 and 140 m respectively.
Figs. 7 and 8 are used to determine the spatial limits of the Monte-Carlo simulations. Equation (10) is used to establish the depth profile in the ground. Using these analytically established parameters as the limits for a Monte-Carlo model will result in a accurate prediction of the response of a gamma-ray spectrometer at heights up to 40 m.

Monte-Carlo results
Monte Carlo simulations for heights up to 40 m with an interval of 5 m have been done. To establish the spatial limits of the geometry of the input model, the methods described in section 3.2 and 3.3 are used. Fig. 9a and b shows the broadened spectra for a selection of these Monte-Carlo simulations for 40 K and the 232 Th-series. From these figures it can be seen that all the peaks present in the spectra decrease with increasing detector height.
The decay of 40 K only emits a single gamma ray with an energy of 1.46 MeV and this is reflected in the shape of the spectrum: a single full energy peak and a Compton continuum. Fig. 10 shows the peak to Compton ratio of the 1.46 MeV (integrated between 1.37 and 1.57 MeV) peak and the integrated Compton continuum between 0.3 and 1.36 MeV. A decrease of 25% in the peak to Compton ratio is observed when Fig. 7. Relative intensity as a function of ground volume for two different heights and ground depths. The volume increase is realized by an increase of the angle in equation (11). moving the detector from 0.80 to 40 m height.
The 232 Th-series consists of a decay chain emitting 346 different gamma energies. The spectrum in Fig. 9b is the sum of all the individual gamma energies and their associated Compton continuum. Fig. 10 shows the relative ratio between the 2.62 MeV peak and the peak at 0.58 MeV. Both peaks result from the decay of 208 Tl. The peak-to-peak ratios for both the analytical and the Monte Carlo models are shown in the figure. The Monte-Carlo simulations result in a peak-to-peak ratio that increases with 35% at a height of 40 m.
It should be noted that the reason for the decrease of the peak-to-Compton ratio of 40 K is the increased number of scattered photons with increasing height. The peak-to-peak ratio in the thorium spectrum is governed by the relative larger attenuation of photons in the low energy peak, with increasing height. These are two significantly different effects.
The intensity of the characteristic gamma energies of 40 K (1.46 MeV) and the 232 Th-series (2.62 MeV), commonly used in windows analysis, are plotted in Fig. 9c and d. The intensities have been extracted from the unbroadened spectra. In the same plot a fit of equation (2) has been added, which has been fitted to this data by using the prefix (Aεγ) as a free parameter and using the associated values for h, μ a , μ g , ρ a , and ρ g . It can be observed that this fit is a good description of the reduction of the peak intensity that result from the MCNP simulation.

Analytical height corrections vs Monte-Carlo simulated spectra
Fig. 9c and d shows that the reduction in peak intensity is described by the exponential integral given in equation (2). To study the spectral shape changes as a function of height, ground based standard spectra at 0.80 m height were scaled to the simulated altitudes. The spectra that result from a Monte-Carlo simulation with the detector at a certain reference height h ref above the ground can be scaled to higher altitudes by using the analytical result given by equation (2). The intensity in each channel is scaled accordingly: where I total is given by equation (2). S ref i is the intensity of the i-th channel of the spectrum at h ref . Each channel corresponds to an energy E i , as defined by the input reference spectrum. Equation (13) contains the ratio of I total and the prefix in equation (2) disappears. Therefore, S scaled i is independent of the detector characteristics and soil properties, but still contains an energy dependent value of μ a (E i ) that is embedded in the exponential integral. Equation (13) can be implemented by using a value μ a (E i ) that is energy dependent as defined by the i-th channel of the reference spectrum. A second approach is implemented by using a fixed value for μ a which corresponds to the dominant energy present in the spectrum. For 40 K the μ a (1.46 MeV) is used and for 232 Th μ a (2.62 MeV) is used. Fig. 11 shows the resulting spectra (in blue) when equation (13) is used to iterate over the channels (i) of a reference spectrum (S ref ) at 0.80 m height to scale the reference spectra to the heights up to 40 m for the two approaches. In red are the scaled spectra that use a fixed μ a to scale the spectra. Both scaled spectra are compared to the spectra resulting from a Monte-Carlo simulation at the target height (in green). Fig. 11a shows the spectra for 40 K at three different heights. It can be seen that the intensity, for all heights, of the main energy peak in the Fig. 8. Radius of the analytical model as a function of height. The lines represent radii that include 65%, 95%, and 99% of the intensity for a gamma photon with energy 2.62 MeV. scaled spectra coincide with the intensity of the simulated spectra. However, the Compton continuum of the scaled spectra underestimates the true intensity of the Compton continuum that results from a Monte-Carlo simulation. This mismatch is larger for the spectra that are scaled with μ a (E i ), compared to the scaling of the whole spectrum with μ a (1.46 MeV). The mismatch becomes larger with increasing height. Fig. 11b shows the scaled spectra that result from the 232 Th-series. The highest peak of both scaling methods coincides for all heights with the Monte-Carlo simulations. Lower energy peaks are either underestimated (for the scaling with μ a (E i )) and overestimated for the scaling with μ a (2.62 MeV). Likewise, as for the 40 K scaled spectra, it is observed that the Compton continuum between the peaks is underestimated, for both scaling methods. Fig. 12a and b show spatial concentration maps for 40 K, and 232 Thseries resulting from the ground-based survey. The radionuclide concentrations (average of the topsoil ±1σ) for the field are 690 ± 93 Bq/kg and 35 ± 10 Bq/kg for 40 K and 232 Th respectively.

Height measurements
In these maps the position where the height measurements have been done is marked with a +. Around this center point, three circles with radii 25, 150 and 300 m have been drawn which approximate the origins of 65%, 95% and 99% of the detected radiation when the detector is measuring at a height of 20 m (Fig. 8). A forest on the west side and a road on the right side of the field determined the spatial measurement boundaries.
The two radionuclides have a different absolute concentration. To compare the spatial variation of the two nuclides with each other, both maps have been plotted with limits that extend to ±50% of the average value. It can be seen that the largest spatial variation is found in the thorium concentration (Fig. 12b). For the potassium concentration, most spatial concentration variation lies within 1σ (Fig. 12a).
The 40 K and 232 Th concentrations that result from the measurement taken at the heights 13, 17, 21, 26, 31 and, 35 m have been plotted in Fig. 12c and d. The point at 0.80 m has been retrieved from the groundbased measurements by summing all spectra taken within a radius of 25 m around the center point. Table 1 shows the radionuclide concentrations of the ground measurement and the average of the height measurements.
In Fig. 12c and d it is observed that the resulting nuclide concentrations from the full spectrum analysis with the ground-based spectra decrease with height. Standard spectra were generated by interpolating Monte-Carlo generated spectra to the measurements height. The green area in these figures represent the average concentration ±1σ of the measurements done in the air. Taking the ±1σ error margins into account, the resulting concentrations can be approximated by a straight line. This straight line is an indication that the radionuclide concentration can be determined independent of the measurement height.
A Z-test to compare the measurements at height versus the ground measurement results in P values of 0.22 and 0.41 for 40 K and 232 Th respectively. These values for P are well above the conventional decision limit of 0.05 and indicate that there is no significant difference in radionuclide concentration between the ground and height measurement.

Discussion
In this study, we presented the footprint and height corrections for gamma-ray spectra associated with the decay of 40 K and 232 Th-series as a function of measurement heights up to 40 m above the ground. It is shown that the reduction in peak intensity of the spectra can be modeled by the exponential integral of the second order, given by equation (2). This exponential integral is used to scale spectra that are simulated for a height of 0.80 m to higher altitudes. These scaled spectra for various heights are compared to spectra that resulted from Monte-Carlo simulations. Moreover, field measurements have been done starting from 0.80 up to 35 m. These measurements have been analysed with and without correcting for height to extract the radionuclide concentrations in the ground.

Part 1: analytical
We presented a new method to calculate the footprint of the radiation in a gamma-ray detector that originates from the ground below by calculating the contribution that falls within a constant attenuation isoline in the ground. The result of this analytical approach expands on the previously published approaches (Duval et al., 1971;Grasty et al., Fig. 10. Peak (integrated between 1.37 and 1.57 MeV) to (0.3-1.36 MeV) Compton ratio resulting from 40 K decay and the 2.62-0.58 MeV peak to peak ratio resulting from 232 Th day. The black line represents the theoretical shift in peak-to-peak ratio based on equation (2) for the 2.62-0.58 MeV peaks. Fig. 11. Scaled spectra compared to Monte-Carlo spectra. a) for 40 K. b) for the 232 Th series. The green lines are spectra resulting from a Monte-Carlo simulation at the indicated height. A Monte-Carlo reference spectrum at 80 cm is scaled to the target height by using equation (13) with an energy dependent μ a (blue line) and a fixed μ a (red line). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) 1979). The expansion comprises a properly defined, finite source volume to be included into the calculation. This finite volume is determined by an isoline that represents the surface from which gamma-rays originate that have an equal intensity contribution to the signal in the detector. It was found that the largest (and deepest) part of the contributing ground volume is located directly below the detector and the depth quickly decreases with increasing angle θ between the detector and ground section of interest. The tails of the isoline contribute only a small fraction of the total radiation intensity, as shown in Fig. 7.
The analytical model only accounts for the reduction of full energy gamma-rays, it does not predict the intensity or shape of the Compton continuum and the 0.511 MeV pair production peak. Therefore, this model does not correctly predict the shape of the spectrum and can only be used to place an upper limit on the size of the area that contributes to the signal in the detector.
The maximum volume that contributes to the signal in the detector is shown in Fig. 7. A plot as shown in Fig. 6 visually represents the origin of the signal. This figure can be used to estimate the footprint (and volume) that contributes to the measured data when the detector is placed at a certain height. This information should be used when designing geophysical surveys using gamma-ray spectrometers, or to establish the geometry of a Monte-Carlo model that predicts the shape and intensity of the gamma-ray spectra.
The altitude of the detector determines the footprint, and thus the flying height determines the spatial accuracy that can be attained. Measurement height should be chosen such that the footprint is smaller than the desired spatial accuracy. The footprint estimation presented in this paper is based on the ground being homogeneous in density and radionuclide content. In real-world situations, where this assumption is not valid, this estimation should be used as a first approximation. Area specific parameters such as the local density can be used to tailor the model.

Part 2: Monte-Carlo simulations
In the second part of this study, the analytical model was used to determine the source-geometry for a Monte-Carlo simulation that would include 99% of the theoretical amount of radiation seen by the detector in a UAV-borne gamma-ray study for altitudes up to 40 m. It was found that a model with a 450 m radius (Fig. 8) and a depth of 75 cm (Fig. 6) is sufficient to model a spectrum that is a good representation of the actual measurement. With current generation of computers, a model running for about 24 h has a statistical precision <5% for the dominant peaks present in the spectrum. The intensities of the 1.46 MeV ( 40 K) and 2.62 MeV ( 232 Th) peaks in the gamma-spectra that result from these Monte-Carlo simulations decrease with detector height, as described by the exponential integral given in equation (2) (Fig. 9c and d).
The Peak to Compton ratio for 40 K reduces by as much as 25% in the height range of 0.80-40 m (Fig. 10). The observation that ground-based spectra that are scaled with μ a (1.46 MeV) agree more with the spectra that result from the Monte-Carlo simulation than the spectra that are scaled using μ a (E i ) is an indication that a large part of the Compton continuum is formed in, or in the vicinity of, the detector. That isif the Comptons were predominantly generated far away from the detector, in the soil, then the absorption would be much better described by absorption coefficients matching their Compton energy.
The number of 1.46 MeV photons, and the accompanying in-thevicinity-formed Compton scattered photons, that end up in the detector scale with equation (2). The extra Comptons that are formed due to the increased layer of air between the source and detector are not accounted for by this equation. The extra Comptons are deduced from the difference between the Monte-Carlo spectra and the μ a (1.46 MeV) scaled spectra. These differences are mostly caused by the extra Comptons that are formed due to this increased layer of air.
The peak-to-peak (0.583 vs 2.62 MeV) ratio for 232 Th day chain increases with 35% in the range of 0.80-40 m (Fig. 10), which indicates that the reduction of the 0.583 MeV peak is sharper than the 2.62 MeV peak. This increase in ratio means that the low energy peaks are attenuated more strongly than the high energy peaks, as predicted by μ a (E i ) in equation (2). However, as can also be seen from Fig. 10, this observed increase in peak-to-peak ratio is smaller than theoretically expected from equation (2). This mismatch means that the underlying Compton continuum must increase with height. This is further substantiated in Fig. 11b in which it can be observed that the Compton continuum in between the peaks behaves similarly to the Compton continuum observed in the 40 K spectra.
The scaling of the ground-based spectra with μ a (2.62 MeV) is a better match for the Compton continuum in between the peaks, but this overestimates the peak intensities of the lower peaks (Fig. 11b). The   Fig. 12. Maps of the agricultural field displaying the interpolated concentration in the soil of: (a) 40 K (average value for the field = 690 ± 93 Bq/kg); and (b) 232 Th (average value for the field = 35 ± 10 Bq/kg) based on data collected in the ground survey. Concentrations at different flying heights (above the location marked with a cross). Both maps have been plotted with limits that extend to ±50% of the average value. (c) 40 K; and (d) 232 Th, obtained using standard spectra at 0.80 m (circles) or interpolated standard spectra at the measurement height (squares). The shaded area represent ±1σ around the average value of the height measurements.

Table 1
Radionuclide concentrations extracted from Fig. 12. Note that these concentrations at 0.80 m is higher than the average concentration of the whole field. This is in agreement with Fig. 12a and b which shows that the point of the height measurement is located above an area with relatively high radionuclide concentrations. overestimation of the peaks in the lower half of these μ a (2.62 MeV) scaled spectra is because the vicinity effect that justifies the use of a μ a for this particular energy cannot be applied if there are multiple initial photons of different energies that make up the spectrum, as is the case for the 232 Th-series. The overestimation of μ a (2.62 MeV) scaled spectra and the underestimation of μ a (E i ) scaled spectrum indicate that two processes are at play that make up the shape of the spectrum. The first process is the reduction in peak intensity of all the peaks present in the 232 Th-series that is described by equation (2). The second process is an increase in Compton continuum that can be seen in the 40 K spectrum. The Compton continuum is the sum of all the Compton continua of the individual peaks present in the 232 Th-series. These two processes describe the shape of the spectrum as a result from the 232 Th-series as a function of height. This discrepancy of Compton continuum between the scaled and simulated spectra ( Fig. 11a and b) implies that for both the windows method as for the full spectrum analysis to extract the nuclide concentrations in the soil, a Monte-Carlo based correction should be applied. However from the same figures it is concluded that small height differences can correctly be predicted by using equation (13) to scale spectra. However, the mismatch increases with increasing difference between the reference spectrum (h ref ) and the target scaling height (h). This leads to the conclusion that Monte-Carlo based results, for a selection of heights, should be used for the analysis of UAV borne gammaray measurements. For the windows method these simulations are used to predict the stripping ratios while for the full spectrum analysis approach the shape of the whole spectra are needed as this method relies on standard spectra. Fig. 12 shows the results of an actual measurement taken on the ground and at several heights. The measurements have been analysed with and without correcting for height to extract the radionuclide concentrations in the ground. From the spatial nuclide distribution maps 10a-b it can already be concluded that the large variation for the 232 Th concentrations could interfere with the concentration determination. When moving to higher altitudes, a larger area is included in the measurement which will contribute to an average radionuclide concentration. Compared to the 232 Th concentration, there is only a small variation in the concentration of 40 K throughout the field. Therefore, it is expected that the concentration prediction for 40 K at height will be more accurate than for 232 Th.

Part 3: height measurements
The 40 K concentrations extracted from spectra taken at various altitudes are seen to decrease if a single 0.80 m calibration is used. However, when extracted using the proper elevated standard spectra, the concentrations come out constantas required. There is a small variation observed in these concentrations at height that is likely to be caused by the sum of the differences between the heights for measurement and simulation, and the spatial inhomogeneity of the underlying area.
A similar result is obtained for the 232 Th data analysed with the height spectra. As with the 40 K concentration the analysis with the ground based spectra decreases with height. Compared to the 40 K, this 232 Th concentration exhibit more variation. One of the factors that influences this result (more so than for 40 K) is the spatial inhomogeneity in 232 Th concentration of the underlying area. The models used to analyse the spectra are based on homogenous concentration distribution. Therefore, it is concluded that the flying height is critical for capturing the true spatial variation in nuclide concentration in an area. Flying at a height where the footprint (Fig. 6) exceeds the desired spatial accuracy will result in the average concentration of the ground that lies within the footprint volume.
The concentrations at 0.80 m, for both 40 K and 232 Th, matches the concentrations that result from the measurements done at the different heights. This means that analysing gamma-ray spectra measured at altitudes between 0 and 40 m can be used to find accurate absolute radionuclide concentrations in the soil.
The radionuclide concentrations shown in Fig. 12c and d are determined by using full spectrum analysis and using interpolated standard spectra. It is concluded that a finite set of simulated spectra (at an interval of 5 m height) is sufficient to approximate standard spectra for heights up to 40 m. With the method used in this research it is hard to distinguish if the variations in the concentration versus height is caused by the standard spectra or by the underlying radionuclide distribution of the field, since both parameters vary with height.

Conclusion
Height corrections for gamma-ray spectra are applied by using the exponential integral. It is shown by this analytical approximation, and also by the Monte-Carlo validation, that the attenuation of gammaradiation as a function of measurement height significantly influences the measurement, with a reduction as much as 50% in intensity at a height of 40 m. This degree of attenuation is only valid for the full energy peaks and the exponential integral does not account for an increased number of scattered photons that end up in the Compton continuum. The shape of this Compton continuum is difficult to predict by an analytical model, however it can be predicted by a Monte-Carlo simulation. This has been shown and verified in this paper for the range of 0-40 m height.
A finite set of Monte-Carlo simulated spectra, at an interval of 5 m, can be used to create standard spectra for gamma-ray measurement in the range of 0-40 m altitude. These standard spectra should be used to determine radionuclide concentrations with full spectrum analysis (which uses these standard spectra directly), or can be used to determine the stripping ratios for the windows method.
It is well known that the footprint of a gamma-ray measurement increases with measurement height and determines the obtainable spatial resolution of an airborne survey. In this work a model that characterizes the origin of radiation as a function of spatial position in the ground is presented to more accurately predict what parts of the ground contribute to the signal in the detector.
It is recommended that the flying height of a UAV-based gamma-ray survey is determined by a footprint that includes between 65 and 95% of the radiation (Fig. 8). This means that the smallest achievable footprint at 10 m flying height lies between 22 and 91 m and between 40 and 140 m at 20 m flying height.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests. The author S. van der Veeke is employed by Medusa Radiometrics and sponsored with an EFRD grant from the European Union. The author declares no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results