Sensitivity Analysis of the DART Model for Forest Mensuration with Airborne Laser Scanning

: Airborne Laser Scanning (ALS) measurements are increasingly vital in forest management and national forest inventories. Despite the growing reliance on ALS data, comparatively little research has examined the sensitivity of ALS measurements to varying survey conditions over commercially important forests. This study investigated: (i) how accurately the Discrete Anisotropic Radiative Transfer (DART) model was able to replicate small-footprint ALS measurements collected over Irish conifer plantations, and (ii) how survey characteristics inﬂuenced the precision of discrete-return metrics. A variance-based global sensitivity analysis demonstrated that discrete-return height distributions were accurately and consistently simulated across 100 forest inventory plots with few perturbations induced by varying acquisition parameters or ground topography. In contrast, discrete return density, canopy cover and the proportion of multiple returns were sensitive to ﬂuctuations in sensor altitude, scanning angle, pulse repetition frequency and pulse duration. Our ﬁndings corroborate previous studies indicating that discrete-return heights are robust to varying acquisition parameters and may be reliable predictors for the indirect retrieval of forest inventory measurements. However, canopy cover and density metrics are only comparable for ALS data collected under similar acquisition conditions, precluding their universal use across different ALS surveys. Our study demonstrates that DART is a robust model for simulating discrete-return measurements over structurally complex forests; however, the replication of foliage morphology, density and orientation are important considerations for radiative transfer simulations using synthetic trees with explicitly deﬁned crown architectures.


Introduction
Forest monitoring programs have been revolutionized by the growing availability and integration of remotely sensed data [1][2][3], with Airborne Laser Scanning (ALS) at the forefront of research in forest management and national forest inventories [4][5][6][7]. Numerous studies have now demonstrated that ALS height and density metrics (calculated at the stand or single-tree scale) are correlated with forest inventory measurements recorded at ground level. These include biometrics such as tree height [8][9][10], stem diameter [11,12], canopy cover and gap fraction [13,14], leaf-area index [15,16], timber volume [17,18], biomass [19,20] and carbon content [21,22]. Consequently, researchers and forest managers have increasingly used ALS measurements for the indirect retrieval of forest biometrics that are impractical to quantify directly through conventional forest mensuration techniques [23][24][25][26][27]. Forest inventory measurements derived from ALS data are often extrapolated to regional or national scales using covariates obtained from optical or radar satellite observations [28][29][30]. This multi-sensor upscaling approach facilitates the wall-to-wall retrieval of forest inventory biometrics through the integration of wide-area satellite observations, high-resolution ALS data and field measurements acquired from a representative sample of forests [31][32][33][34].
Today, most ALS sensors used in operational forestry are small-footprint (<1 m) pulsed laser systems capable of recording a finite number of discrete returns (typically <5) per laser pulse [35]. Since discrete-return instruments decompose the backscatter signal into a limited number of targets, they provide vestigial measurements of forest structure and are often insensitive to stem and understory components occluded by the forest canopy [36,37]. In contrast, waveform systems digitize the entire backscatter signal at high frequency, thereby recording more detailed information on forest structures residing within the laser footprint [38,39]. While waveform systems have become increasingly prevalent in recent years, ALS surveys undertaken for operational forestry primarily employ discrete-return metrics obtained through the Gaussian decomposition of ALS waveforms [40,41]. This may be due to the additional cost and complexity of processing waveform data [42], or in some cases, due to the marginal gains provided by full-waveform analysis [43][44][45][46].
Forest biometrics are derived from discrete-return measurements using two common methods, namely the area-based approach (ABA) [6,47,48] and the single-tree approach [49][50][51]. In the ABA, ALS measurements are spatially gridded and converted into a range of structural measures (height percentiles and density statistics) that quantify the vertical and horizontal distribution of laser returns within the forest canopy [52][53][54][55]. Discrete-return statistics are then regressed against forest inventory measurements recorded at ground level in fixed-area plots [56][57][58][59]. While the ABA facilitates the indirect retrieval of forest inventory measurements with relatively inexpensive low-density ALS data [60][61][62], the method explicitly assumes that variations in the distribution of laser returns are directly related to forest structural properties (e.g., tree height, crown morphology, foliage density), when in fact, discrete-return distributions are also influenced by variations in sensor configuration, viewing geometry and ground topography [63][64][65]. Such variations can impact the derivation of forest biometrics, causing errors of ±10% for timber volume retrieval when predictive models are not reconstructed or refitted for different ALS surveys [66,67].
Previous studies have examined the influence of acquisition parameters on the precision of ALS measurements through repeat surveys undertaken over individual research sites [68][69][70][71]. Several such investigations have demonstrated that discrete-return metrics are influenced significantly by variations in platform altitude, scanning angle and pulse repetition frequency [66,67,72,73]. Higher pulse repetition frequencies are associated with a reduction in pulse energy that decreases pulse penetration [74,75] and the proportion of multiple canopy returns [66,73]. Conversely, lower pulse frequencies and increased platform altitude are associated with a reduction in pulse density that results in canopy height underestimation as: (i) fewer returns are recorded directly at treetops [63,76,77], and (ii) fewer ground returns are acquired for the interpolation of digital terrain models with accurate microtopography [73,78,79]. The influence of pulse density has also been tested empirically with different sensors [74,80,81] and through the artificial decimation of discrete returns [61,82,83]. Several studies have established that low-density ALS data (0.5-2 pulses m 2 ) are sufficient to obtain reliable height and density metrics for the retrieval of most forest inventory parameters [60,62]. However, forest biometrics such as canopy cover, gap fraction and leaf-area index are often measured less precisely due to their greater sensitivity to variations in discrete-return density [61,84] and laser scanning angle [85,86].
While empirical studies can provide valuable insights into the impact of survey characteristics on the precision of ALS measurements, repeat surveys do not offer ideal experimental conditions due to: (i) the anisotropy of ALS backscatter arising from variations in sensor position and viewing geometry, and (ii) sensor-specific characteristics (e.g., triggering mechanism, scanning pattern, laser beam profile) that complicate the task of isolating the effects of individual acquisition parameters. An alternative approach involves the testing of acquisition parameters through radiative transfer simulations where ALS surveys are replicated in controlled experiments over computer-generated forests [65,87,88]. Radiative transfer simulations provide a cost-effective alternative to repeat ALS surveys, while allowing a greater variety of acquisition parameters to be tested systematically across a broader range of forests than can be achieved through direct measurements [89]. Moreover, the information gleaned through the physical modeling of ALS data can enhance our understanding of how laser pulses interact with forest canopies [90][91][92][93], enabling the retrieval of forest biometrics through radiative transfer inversion techniques without requiring exhaustive field measurements needed for empirically derived predictive models [14,[94][95][96][97]. Despite these advantages, the task of replicating forest canopies in radiative transfer models remains an active area of research due to challenges associated with: (i) reproducing individual tree locations, (ii) generating synthetic trees with realistic structure and allometric properties, and (iii) selecting crown archetypes that adequately represent forest canopy structure and reflectance characteristics [98][99][100][101][102].
While several simulation studies have examined the impact of sensor and survey design on large-footprint (≥3 m) waveform LiDAR (Light Detection and Ranging) acquisitions [93,94,97,[103][104][105], fewer investigations have examined small-footprint discrete-return measurements due to the computational demands of simulating multi-pulse ALS acquisitions over complex landscapes [65,96,[106][107][108]. In previous studies, the influence of survey characteristics on ALS measurements has been tested through local sensitivity analyses, where only one acquisition parameter was varied at a time while the remaining parameters were held constant at nominal values [65,87,88,109,110]. Unfortunately, one-at-a-time methods are only appropriate for additive or linear models since they: (i) inadequately sample the parameter feature-space, and (ii) are unable to account for non-linear interactions between the simulation parameters [111][112][113]. In addition, forest canopies are often represented using turbid layers or geometric primitives (conical and ellipsoidal objects comprised of turbid media) as proxies for individual tree crowns [94,[106][107][108][114][115][116][117]. While the use of abstract forest canopies can produce accurate simulations in certain forests [118,119], recent research demonstrates that simple crown archetypes can introduce large inaccuracies when simulating reflectance over structurally diverse forests due to the inadequate replication of twig and foliage morphology, density and orientation [98,100,[120][121][122]. These findings indicate that more detailed three-dimensional representations of individual trees, e.g., [102,123,124] are required for realistic radiative transfer simulations over heterogeneous forests.
In the present study, we attempt to avoid a perfunctory sensitivity analysis [111,113] of the Discrete Anisotropic Radiative Transfer (DART) model by: (i) simulating ALS acquisitions over synthetic forests with explicitly defined crown architectures, and (ii) performing a variance-based sensitivity analysis in which all survey parameters are modulated systematically to identify factors influencing the accuracy and precision (repeatability) of discrete-return metrics [125,126]. Our investigation focuses exclusively on the simulation of small-footprint discrete-return ALS data due to its widespread use in operational forestry and its growing availability following the deployment of discrete-return sensors onboard unmanned aerial vehicles [127,128]. The aims of this study are to: (i) investigate whether a radiative transfer model can accurately replicate discrete-return measurements collected over commercial forests in the Republic of Ireland, and (ii) explore how survey and model parameters influence the precision of discrete-return metrics simulated over structurally diverse forests.

The Discrete Anisotropic Radiative Transfer (DART) Model
Radiative transfer in forested environments can be simulated with semi-empirical analytical models [129][130][131], Geometric Optical (GO) models [132][133][134] or radiative transfer models where the propagation of electromagnetic radiation is replicated through physically based mechanisms that govern photon interception, scattering and transmittance [88,135]. Analytical radiative transfer models commonly employ closed-form mathematical functions to replicate canopy reflectance measurements obtained from field or laboratory experiments, consequently their application is often limited by: (i) simplified representations of radiative transfer mechanisms, (ii) a lack of generality in their empirical parameterization, and (iii) an inability to simulate heterogeneous forests with complex ground topography [89,136,137]. A more detailed, flexible and transferable approach to radiative transfer modeling can be achieved through ray-tracing experiments where photon trajectories are simulated through Successive Orders of Scattering [138], radiosity methods [139,140], discrete ordinates [141,142] or Monte Carlo techniques [89]. Monte Carlo ray-tracing methods provide a tractable solution to radiative transfer modeling in forest canopies where photon scattering is stochastic, requiring that photon trajectories be repeatedly sampled in order to characterize bi-directional reflectance across different sensor viewing geometries and operating wavelengths. Monte Carlo methods have gained in popularity and complexity due to the recent proliferation of computer resources, with several models designed to simulate radiative transfer in one-dimensional canopies (e.g., 4SAIL [143], 5-Scale [144,145]) and three-dimensional forests (e.g., FLIGHT [87,146], Raytran [114], librat [147,148], FLiES [149]).
Three-dimensional radiative transfer models are required when path length is explicitly simulated for ALS acquisitions or radar interferometry [150]. DART is a three-dimensional model designed to simulate the propagation of electromagnetic radiation through heterogeneous terrestrial landscapes while accounting for variations in solar irradiance, atmospheric composition and sensor viewing geometry [137,[151][152][153]. DART simulations of vegetation reflectance in the visible and near-infrared wavelengths have been used for the inversion of vegetation biometrics [154][155][156][157][158][159][160] and for reproducing imagery acquired from satellite and airborne sensors [161][162][163]. DART simulations have also been validated against field measurements [164] and compared against several independent reflectance models in the Radiation transfer Model Intercomparison (RAMI) experiments [165][166][167][168][169]. Although DART has been tested extensively in the passive optical domain, its ability to simulate small-footprint ALS data has received less attention, e.g., [170,171], with previous work focusing on terrestrial laser scanning [172] and large-footprint waveform LiDAR [173][174][175]. However, recent improvements to DART have facilitated the efficient simulation of multi-pulse acquisitions [137,176] with the ability to export ALS waveforms into the Sorted Pulse Data (SPD) format [177].
DART simulates the propagation of laser pulses in terrestrial landscapes using a forward Monte Carlo ray-tracing framework where pulse trajectories are tracked from the sensor using physically based mechanisms that govern photon interception, scattering and thermal emission [137,178]. To facilitate rapid and computationally efficient simulations, DART employs an Iterative Uniform Squared Discretization (IUSD) scheme to constrain pulse tracking to a finite number of pre-defined discrete scattering directions (Ω D ) that evenly sample the bi-directional reflectance distribution [179,180]. This pulse tracking framework greatly accelerates DART simulations as probable pulse trajectories can be determined from pre-computed scattering matrices using the box method [173]. However, this simplified scattering framework is susceptible to discretization error [181,182] if reflectance is anisotropic and the number of pre-defined angular directions are too few to accurately characterize the bi-directional reflectance distribution [178]. Such errors are negated by specifying many pre-defined scattering directions to produce a finer sampling of bi-directional reflectance at each scattering event [180].
DART can simulate terrestrial landscapes at varying levels of complexity depending on: (i) the level of detail required when replicating remote sensing measurements, and (ii) the computational resources available for performing ray-tracing simulations. Terrestrial landscapes are simulated in a voxel matrix where forest canopy components (stems, branches, foliage) are represented as triangular surfaces, turbid media or a combination of both. Turbid media voxels provide a simple representation of plant biomass, comprising of small planar surfaces with uniform density, angular distribution and reflectance properties. Conversely, for complex scenes where the position and geometry of forest canopy components are explicitly defined, DART voxels are comprised of triangles and parallelograms with discrete reflectance and transmittance characteristics [137,178]. In this latter mode of operation, DART can replicate structurally complex forests through the direct importation of 3D trees generated in external modeling software [183,184] or trees digitized from terrestrial laser scanning measurements [124,185]. Recent work by Roberts et al. [186] and Janoutová et al. [187] suggests that forest biometrics may be reliably inverted using simplified trees with foliage represented by turbid media. However, in the present study, we seek to minimize geometric assumptions associated with simplified crowns, therefore ALS acquisitions are simulated using trees with explicitly defined crown geometries and reflectance characteristics. We anticipate that trees with detailed crown architectures will enable more accurate ALS simulations at the scale of individual forest inventory plots, thereby permitting a robust sensitivity analysis of the DART model to variations in sensor and survey characteristics.

Reference Data
To calibrate our DART simulations, an ALS survey was acquired over 27,000 hectares of Sitka spruce (Picea sitchensis) and Lodgepole pine (Pinus contorta) plantations in the west of Ireland ( Figure 1). These forests were selected due to their structural heterogeneity, with some compartments consisting of unthinned juvenile stands, while other compartments contained mature first-rotation crops having undergone several thinning treatments. The ALS survey was undertaken in February 2015 with a RIEGL LMS-Q680i sensor onboard a Pilatus PC6 aircraft flying at a nominal altitude of 550 m and a mean velocity of 204 ± 5 kph. 116 flight lines were flown with a 40% overlap, a pulse repetition frequency of 360 kHz and a maximum scanning angle of 30 • from nadir (60 • swath), producing a mean point density of 29 ± 10 returns m 2 . The ALS data was processed with SPDLib [188] for the purpose of ground classification, the calculation of normalized heights and the production of area-based discrete-return metrics at 20 m resolution (cell size = 400 m 2 ) throughout the study region. Ground classification followed a two-stage procedure where an initial ground surface was derived through a progressive morphology filter [189] and ground returns were then refined using a multi-scale curvature algorithm [190]. Normalized vegetation heights were calculated from ground surface elevation interpolated using the natural neighbor algorithm [191].  To accompany the ALS survey, Coillte Teoranta (Irish State Forestry Board) surveyors undertook a forest inventory campaign, collecting 170 fixed-area (400 m 2 ) plots distributed throughout the study region ( Figure 1). Plot centroids were geo-referenced with a differential GNSS receiver and the diameter of each tree was measured at breast height with digital calipers. Within each fixed-area plot, five tree height measurements were collected with a Vertex IV clinometer (www.haglofcg.com) for trees corresponding to the minimum, mean, maximum, 25th and 75th percentiles of the stem diameter distribution for the dominant tree species. Height measurements were also collected for non-dominant tree species (e.g., Larix kaempferi, Pseudotsuga menziesii, Betula pendula), producing DBH and height measurements for a total of 1232 trees. The forest inventory revealed structural heterogeneity in the plantations, with stocking densities of 525-2375 trees ha −1 , basal area varying between 6 and 82 m 2 ha −1 and tree heights ranging from 2.5-30.9 m (Table 1).

Forest Plot Reconstruction
Allometric relationships were derived from the single-tree measurements to facilitate the creation of 700 computer-generated trees using the OnyxTree modeling software [192]. A diverse range of trees were generated with heights ranging from 3 to 30 m, DBH values of 3-73 cm and crown diameters of 1.2-8.6 m. Each simulated tree was comprised of a single stem connected to a series of boughs, branches, twigs and needles. However, the position, length and orientation of each component was randomly varied to ensure that each tree was structurally unique. Needles were simulated as rectangular planar surfaces ( Figure 2c) with a length of 30 mm ± 15% and inclined at 75-90 • from each twig to replicate the erectophile needle orientation prevalent in Sitka spruce [193]. Stem profiles were modeled using a taper function developed for Sitka spruce plantations in Northern Britain: where z represents the stem height as a fraction of top height, dbh is the tree diameter at breast height and d denotes the estimated tree diameter at height z [194]. To accelerate the DART simulations, tree components were organized into three structural groups (stems, branches and needles) and Lambertian reflectance properties were assigned to each group based on field spectroradiometer measurements collected over the 350-2500 nm spectral range. 100 fixed-area (400 m 2 ) plots from the forest inventory campaign were selected for reconstruction based on a maximin Latin hypercube sampling strategy which ensured that the reconstructed plots were sampled throughout the multidimensional distribution of the forest inventory attributes [195][196][197]. Individual tree positions and biophysical attributes (e.g., height, DBH, crown diameter) were then estimated through Monte Carlo simulations that sought to minimize differences between the simulated and observed plot biometrics in terms of tree density, basal area and the minimum, mean and maximum values for DBH and tree height. This statistical approach was necessary because: (i) individual tree locations were not recorded in the forest inventory, and (ii) tree detection methods (e.g., watershed segmentation, local maxima and region-growing [198,199], point-cloud segmentation [200,201] and manual digitizing) were unsuccessful in delineating individual trees with sufficient accuracy due to the presence of overlapping crowns and suppressed trees. The synthetic trees were positioned at regular intervals along parallel rows to replicate the planting patterns and stocking densities found in the Slieve Aughty plantations (Figure 2a,b). A random offset was applied to each tree position to prevent perfectly straight planting rows, and realistic ground topography was incorporated by importing digital terrain models from the reference ALS survey. To replicate edge-effects caused by nearby trees protruding into the forest plot, each fixed-area plot was positioned within a 900 m 2 DART scene with additional trees located around the plot perimeter. Following each DART simulation, the ALS data was subset to the central 400 m 2 plot by discarding laser returns recorded outside of the plot boundary.

DART Simulation Parameters
DART (version 5.7.3) simulations were performed in multi-pulse mode using finite scenes where pulses were discarded once they had exited the simulation boundary [176]. To mimic acquisition parameters used in the reference ALS survey, we restricted laser scanning angles to ≤30 • from nadir (60 • swath) and we maintained a constant laser operating wavelength of 1550 nm to replicate the Erbium-doped fiber laser onboard the RIEGL LMS-Q680i sensor. Following [75], laser beam divergence was modulated between 0.3-0.8 mrad to investigate the impact of spreading pulse energy across variable-sized (0.12-1.2 m) laser footprints. To replicate the work of Yu et al. [77] and Takahashi et al. [78], platform altitude was varied from 400-1500 m. Pulse repetition frequency (PRF) was adjusted from 100-400 kHz, while the scan rate was modulated between 50-200 Hz. These parameter ranges were consistent with the RIEGL LMS-Q680i sensor specifications [202], producing pulse densities (8-30 pulses m 2 ) consistently higher than for ALS acquisitions in Scandinavia where sparse data (0.5-2 pulses m 2 ) are often used for area-based forest inventory [66,67].
Several of the DART model parameters (Table 2) were modulated to investigate their influence on discrete-return metrics. The number of pre-defined discrete scattering directions (Ω D ) was varied from 60-200, while the voxel resolution was adjusted between 0.1 and 1.0 m. We expected these two simulation parameters to have a negligible impact on the discrete-return metrics due to the avoidance of turbid media [203]; however we chose to vary these parameters as a sanity check. Each simulated pulse was tracked until its energy fell below a propagation threshold of 5-15% of mean pulse energy, while the discretization interval for digitizing each received waveform (acquisition rate) was varied between 0.5 and 3 ns. The minimum number of simulated photons fired into each DART cell (sub-center photons) was adjusted between 1 and 15 to examine the impact of increasing the number of Monte Carlo iterations performed for simulating photon scattering trajectories within each voxel. The scattering order, controlling the maximum number of scattering events that each simulated photon can undergo, was also modulated between 5 and 15. Laser pulses were fired along parallel scan lines to replicate the scanning pattern produced by the rotating polygon mirror onboard the RIEGL LMS-Q680i sensor [204]. Equations for generating realistic pulse distributions were obtained from Wehr and Lohr [205] and Baltsavias [206]. The along-track pulse spacing (dx along ) was calculated following: where v is the sensor flying velocity (m/s) and f sc is the scan rate measured in hertz or horizontal scan lines per second. The number of pulses per scan line (N pulses ) was then determined: where PRF denotes the pulse repetition frequency in hertz. The swath width was subsequently calculated as: where h denotes the sensor altitude in meters and θ represents the maximum scanning half-angle in radians. Finally, the across-track pulse spacing (dx across ) was determined: thereby enabling pulses to be positioned at regular intervals (i.e., consistent across-track spacing) across each scan line. Based on the empirical observations of Armston [207], peak pulse power (A peak ) was varied as a function of pulse repetition frequency (kHz): and the temporal profile of each laser pulse was simulated using a centered super-Gaussian distribution: where t denotes the elapsed time in ns and σ t represents the time at which pulse intensity reaches 50% of peak amplitude ( Figure 3). We chose to vary the pulse duration at full-width half-maximum between 3.5 and 5.5 ns to investigate the impact of longer pulse durations on the accuracy and precision of the ALS metrics. These parameter bounds were consistent with experiments of Armston [207] in which 10,000 transmitted pulses were recorded from the RIEGL LMS-Q680i sensor with a mean full-width half-maximum of 4.55 ± 0.016 ns. Similarly, Jalobeanu and Gonçalves [204] observed Gaussian pulses with a mean half-pulse duration of 4.3 ns. The Gaussian shape and duration of the transmitted pulse is consistent with prior simulation studies using the RIEGL LMS-Q680i sensor [204,207].
The distribution of laser energy in each pulse footprint was simulated using a two-dimensional super-Gaussian distribution: whereP 1,β denotes the ratio of pulse energy to the central maximum, β is the angular offset from the pulse direction and σ β is the two-dimensional standard deviation of laser energy [203]. The dispersal of laser energy in each pulse footprint was varied by adjusting σ β between 0.2 and 0.5 ( Figure 4). The use of a Gaussian pulse footprint was consistent with the work of Armston [207]; however we could not validate this assumption due to a lack of published research documenting the laser energy profile employed by the RIEGL LMS-Q680i sensor. Pulse intensity (%)

Discrete-Return Metrics
The simulated waveforms from DART were stored in the Sorted Pulse Data (SPD) format [177] and Gaussian waveform decomposition was undertaken in SPDLib [188] with an algorithm designed to replicate the results of RiAnalyze [208]. Normalized discrete-return heights were then calculated by comparing return elevations to the digital terrain surface imported into each DART simulation. Discrete returns were classified as: (i) ground returns if they were within ±1 m of the terrain, or (ii) vegetation if returns were ≥2 m above ground. We also tested several automated ground detection methods in SPDLib (e.g., progressive morphology filter [189], multi-scale curvature algorithm [190], parameter-free filtering algorithm [209]); however a simple threshold proved to be superior to automated methods since we could exploit prior knowledge of the ground surface imported into each DART simulation.
Area-based metrics (Table 3) were calculated at the individual plot scale (400 m 2 ) using custom Python scripts that were multi-threaded with GNU Parallel [210]. Discrete-return metrics included height percentiles often used for the indirect retrieval of forest inventory parameters [52,53]. Canopy cover and gap fraction indices quantified the degree of canopy openness [16], while the relative abundance of single, multiple and ground returns provided a measure of laser penetration rates [211]. Gap fraction was calculated using all return classes (i.e., first, intermediate and last echoes) as this formulation was similar to the All-echo cover index that had proven to be successful in detecting small gaps in closed-canopy forests [212,213]. Canopy cover was calculated using a grid-based binary classification method [13] as the percentage of pixels (at 0.5-meter resolution) containing vegetation returns above 2 m. Thus, canopy cover as defined here represents the proportion of the forest floor occluded by the vertical projection of tree crowns (vertical canopy cover) rather than the proportion of sky obscured by vegetation (angular canopy closure) [214,215]. This method was chosen as it provides an unbiased estimate of vertical canopy cover while being sensitive to variations in discrete-return density and laser scanning angle [13].
Voxel-based canopy volume metrics [216] were also calculated to investigate the vertical and horizontal variability of discrete returns generated by the DART simulations. Voxels (1 m 3 ) above the 65th height percentile (h 65 ) and containing laser returns were classified as Euphotic to represent the overstory where the bulk of incoming radiation was intercepted. Conversely, voxels below h 65 and containing laser returns were defined as Oligophotic to represent the partially obscured understory. The Euphotic and Oligophotic voxels together represented the filled canopy volume, whereas voxels lacking laser returns were classified as Open Gaps or Closed Gaps depending on their relative positions above or below the forest canopy [217,218]. The raw counts of voxel class membership were converted into normalized percentages to facilitate comparisons across all DART simulations. Table 3. Discrete-return metrics calculated for each DART simulation.
h 25 -h max Height percentiles of vegetation returns above 2 m. [48] h mean Mean height of vegetation returns above 2 m.
Multiple return (%) Percentage of laser pulses with multiple returns.
Return density Mean number of discrete returns per m 2 .
Single return (%) Percentage of laser pulses with single returns.

Sensitivity Analysis
Following guidelines set out by Saltelli and Annoni [111] and Saltelli et al. [113], our experiment was divided into two stages: (i) an uncertainty analysis was undertaken to determine the accuracy and precision of discrete-return measurements generated by the DART model, (ii) a sensitivity analysis was conducted to identify the acquisition parameters that had caused greatest variability in the simulated discrete-return metrics.
To achieve the first objective, 500 DART simulations were performed over each of the reconstructed forest plots (50,000 simulations in total) using a sequence of acquisition parameters generated through a maximin Latin hypercube sampling of the parameter bounds identified in Table 2. This sampling approach ensured that: (i) the simulation parameters were sampled uniformly across the parameter feature-space, and (ii) the same acquisition parameters were repeatedly tested over each of the 100 forest inventory plots. We evaluated the precision of the simulated ALS metrics using a reliability ratio: whereσ 2 u represents the between-plot variance andσ 2 w denotes the within-plot variance of each ALS metric [82,219]. A reliability ratio close to 1 indicated that DART simulations had produced precise ALS metrics, whereas values closer to zero demonstrated that changes to the acquisition parameters had produced ALS metrics with a high degree of within-plot variance, thus making these ALS metrics unreliable predictors for the derivation of forest inventory parameters. For comparison purposes, we also calculated the Fisher ratio: to examine the ratio of between-plot to within-plot variance for each ALS metric [220]. Two-sample Kolmogorov-Smirnov and Anderson-Darling tests [221][222][223] were then performed to investigate whether the simulated height percentiles were statistically similar those recorded within the reference ALS survey.
For the second objective, a sequence of 20,000 acquisition parameters were generated to perform a variance-based global sensitivity analysis in SALib [224] using the Sobol method [225][226][227].
The acquisition parameters were tested sequentially over one of the simulated forest plots (illustrated in Figure 2), and the unconditional variance of each ALS metric V(y) was decomposed following: (11) where V i denotes the first-order contribution of the ith model parameter, V ij represents the pairwise interaction of the ith and jth parameters and n is the number of acquisition parameters. First, second, and total order of sensitivity indices were then calculated: where V ∼i denotes the contribution from all simulation parameters excluding the ith parameter. First-order sensitivity indices quantified the output variance directly attributable to parameter i without considering parameter interactions. Conversely, total order sensitivity indices quantified the variance attributable to parameter i when including interactions with other simulation parameters [228,229]. Additive simulation outputs lacking parameter interactions were expected to have identical S i and S T i indices with the sum of total order indices (∑ S T i ) converging at unity. In contrast, strong interactions between the simulation parameters were expected to produce divergent sensitivity indices and ∑ S T i values exceeding unity [111].

Uncertainty Analysis
100 fixed-area (400 m 2 ) plots were reconstructed with biometrics similar to those recorded in the forest inventory campaign ( Figure 5). Tree density was reproduced with a root-mean-square error (RMSE) of 67 trees ha −1 (5.6%), while basal area, mean DBH and mean tree height were replicated with errors of 1.7 m 2 ha −1 (5.4%), 1.4 cm (7.3%) and 0.67 m (4.7%) respectively. These results demonstrated a good correspondence between the field measurements and the biometrics of simulated plots, consequently the DART simulations were representative of stands measured during the forest inventory.
The DART simulations produced discrete-return height percentiles with limited within-plot variance as demonstrated by high Fisher ratios and reliability ratios >0.97 (Table 4). Strong reliability ratios (>0.99) were produced for all height percentiles above h 50 , demonstrating that discrete-return height distributions in the upper forest canopy were highly stable and insensitive to variations in sensor configuration, ground topography or DART model parameterization. Height percentiles in the lower forest canopy were marginally less precise, with reliability ratios of 0.99 and 0.97 for h mean and h 25 respectively.
Weak reliability ratios were recorded for discrete-return density (0.108), canopy cover (0.675), the proportion of multiple canopy returns (0.559) and the proportion of Euphotic (0.588) and Oligophotic (0.367) voxels. The high proportion of within-plot variance for these ALS metrics demonstrated that changes to the acquisition parameters had significantly influenced: (i) the density and horizontal distribution of discrete returns, and (i) the extent of laser pulse penetration through the simulated forest canopies. Despite this, gap fraction and the proportion of voxels classified as open and closed gaps produced comparatively strong reliability ratios of 0.932, 0.896 and 0.913 respectively, indicating that gap fraction indices were relatively stable despite perturbations to the acquisition parameters.  The Anderson-Darling and Kolmogorov-Smirnov tests demonstrated that the simulated discrete-return height percentiles were consistent (at the 99% confidence interval) with those recorded in the reference survey (Table 5). Encouragingly, most DART simulations produced height errors of ±1.5 m (Figure 6a), with a mean absolute error of 0.58 ± 0.46 m (3.3 ± 2.6%) for h max . The mean error for h max across all simulations was −0.31 m, indicating that the top of canopy height was often underestimated and therefore negatively biased. However, discrete-return height percentiles within the upper forest canopy (h 60 -h 99 ) were notably more accurate and precise than those simulated in the lower forest canopy as median height and h 25 were underestimated by −0.8 m (6.3%) and −1.0 m (9.7%) respectively.  The canopy cover and laser penetration metrics generated by DART proved to be highly accurate when compared against the reference survey (Figure 6b). Canopy cover was simulated with a mean absolute percentage error (MAPE) of 3.92%, while the voxel-based canopy volume metrics were simulated with mean errors below 3.85%. However, our DART simulations consistently overestimated the proportion of multiple returns (MAPE = 5.1%), while the number of ground returns was overestimated by 26.1%.

Sensitivity Analysis
The sensitivity analysis revealed that laser scanning angle, pulse repetition frequency and sensor altitude had exerted a significant influence across all of the discrete-return metrics (Figure 7). These three acquisition parameters accounted for 75-83% of the total recorded variance in the discrete-return height percentiles (Figure 7i-l). Scan angle exerted an important influence on discrete return density, canopy cover and discrete return heights (S T i = 0.40-0.46). However, pulse frequency also significantly influenced discrete return density and height distributions throughout the forest canopy (S T i = 0.33-0.36). The voxel-based canopy volume metrics were strongly influenced by scan angle (S T i = 0.36-0.44) and PRF (S T i = 0.30-0.35), highlighting the importance of these two simulation parameters in governing both the horizontal and vertical distribution of the discrete returns.
S c a n a n g l e P R S c a n a n g l e P R F A l t i t u d e P u l s e d u r a t i o n 0.0 0.5 S c a n a n g l e P R F A l t i t u d e P u l s e d u r a t i o n 0.0 0.5 S c a n a n g l e Interestingly, the DART model parameters (propagation threshold, scattering directions, scattering order, sub-center photons and voxel resolution) exerted a negligible influence on the simulated ALS metrics. First and total order sensitivity indices for these simulation parameters were consistently below 0.05 and statistically indistinguishable from zero due to large standard errors. These results demonstrated that DART parameters could be fixed at constant values to reduce the parameterization of future simulations.
We found little evidence of strong interactions between the simulation parameters, as shown by: (i) the homogeneity of total and first-order sensitivity indices (Figure 7), and (ii) the lack of strong second-order sensitivity indices. For instance, sensor altitude, PRF and sensor altitude exerted an important influence on h max (Figure 7l); however we did not observe strong second-order interactions between scan angle and PRF (S ij = 0.072), scan angle and sensor altitude (S ij = 0.05) or between PRF and sensor altitude (S ij = 0.009). The lack of second-order interactions indicated that the discrete-return metrics were the result of additive or multiplicative processes that could be approximated through univariate functions. However, response curves ( Figure 8) revealed a variety of non-linear and non-monotonic relationships, indicating that the discrete-return metrics could not be readily approximated through a linear combination of the simulation parameters. Wider scan angles were associated with a reduction in the discrete-return height distributions ( Figure 8) as a larger proportion of pulses interacted with crown edges rather than tree apices. This downwards shift in the discrete-return heights was most pronounced at wider scanning angles, with height percentiles declining by more than 0.5 m as the maximum scanning angle increased from 20 • to 30 • . We observed a similar trend for PRF, with lower pulse frequencies associated with a monotonic reduction in the ALS height percentiles as fewer returns were recorded at treetops. In contrast, rising sensor altitudes were associated with an upward shift in the discrete-return height distributions, with h max rising by 0.5 m as the platform altitude increased from 400 to 1500 m. For the ALS density and canopy cover metrics, we recorded a decline in the number of discrete returns as PRF was reduced and scan angle and sensor altitude increased. Interestingly, canopy cover decreased by 5% as the scan angle increased from 15 • to 30 • . This marked decline in canopy cover at wider scan angles was an unexpected result caused primarily by the reduced discrete-return density at wider scan angles.
As noted in Section 4.1, the largest errors encountered in our DART simulations were associated with the proportion of ground and multiple canopy returns. Our sensitivity analysis revealed that pulse duration, scan angle and pulse repetition frequency were the main sources of variance for these two metrics (Figure 7b,c). In both instances, the duration of transmitted pulses proved to be the dominant simulation parameter, with longer pulses triggering a monotonic increase in the number of ground returns and a decline in the proportion of multiple canopy returns ( Figure 8). Further experiments revealed that longer pulse durations had caused a decline in the proportion of multiple returns due to the integration of backscatter signals from canopy and ground structures (Figure 9). These results demonstrated that Gaussian waveform decomposition was unable to resolve echoes from near-ground structures as half-pulse duration increased beyond 5 ns.

Discussion
In this study, a variance-based global sensitivity analysis was undertaken to: (i) investigate how accurately the DART model was able to reproduce ALS measurements acquired over commercial forests, and (ii) examine the impact of survey parameters on the precision of discrete-return metrics commonly used for the derivation of forest inventory measurements. While previous simulation studies had examined the influence of survey characteristics using local sensitivity analyses and/or abstract forest canopies, we attempted to undertake a more rigorous sensitivity analysis of the DART model using: (i) three-dimensional synthetic trees with explicitly defined crown geometries and reflectance properties, and (ii) a variance-based global sensitivity analysis that accounted for non-linear interactions between the ALS survey parameters.

Simulation Accuracy
The DART simulations successfully replicated discrete-return ALS measurements collected with the RIEGL LMS-Q680i sensor over commercial forests in the Republic of Ireland. Despite considerable modifications to the acquisition parameters, DART produced discrete-return height distributions that were accurate and precise across 100 forest inventory plots, with a mean absolute error of 0.58 ± 0.46 m (3.3 ± 2.6%) for maximum canopy height. Canopy cover and canopy volume metrics were replicated with mean absolute percentage errors below 4%, demonstrating the horizontal and vertical accuracy of discrete-return distributions produced by the DART model.
The canopy height errors generated by our simulations were consistent with those reported by Suárez et al. [230] and Donoghue and Watt [231] for plantation forests in Scotland and Northern England. Hence, our simulations produced canopy height errors that were comparable in magnitude to direct comparisons between ALS and ground-based canopy height measurements in British conifer plantations. However, we acknowledge that our experiments did not take into consideration common biophysical features of plantation forests such as the presence of forked or wolf trees, nor did we account for height discrepancies induced by tree lean on steep terrain [232,233]. We also acknowledge that the forests used in our simulations were simplified representations of the Slieve Aughty plantations since they were parameterized primarily for Sitka spruce monocultures and did not account for nutrient deficient or waterlogged soils [234], broadleaf encroachment [235], silvicultural treatments (thinning and pruning regimes) and stochastic events such as windblow [236] which regularly impact forest structure and alter the reflectance signal recorded by ALS instruments.
Our simulations consistently overestimated the proportion of multiple canopy returns (MAPE = 5.1%), while the proportion of ground returns was overestimated by 26.1%. These large inaccuracies revealed that the simulated forest canopies were too porous, causing an insufficient number of pulses to be wholly reflected by the forest canopy in order to occlude the ground surface. The lack of ground vegetation in our simulations may have contributed towards more penetrable forest canopies; however it is unlikely that this was a dominant factor given the paucity of understory vegetation in conifer plantations throughout Britain and Ireland. It is also feasible that the elevated number of ground returns recorded in our simulations was symptomatic of ground misclassification during the processing of the reference ALS dataset. However, our experiments suggest that the elevated number of ground and multiple returns were primarily caused by structural inaccuracies in the simulated trees due to the inadequate representation of foliage. This finding is supported by previous studies demonstrating the importance of needle morphology, density and orientation in governing multiple scattering in coniferous forests [237,238]. In particular, Rochdi et al. [239] recorded high levels of multiple scattering when performing radiative transfer experiments with shoots comprised of realistic needles, whereas shoots consisting of abstract needles produced little multiple scattering due to shoot impermeability. Our study therefore highlights the complexities of: (i) generating synthetic trees that are biologically realistic, and (ii) parameterizing tree modeling software such as OnyxTree, Arbaro [240] or AMAPstudio [184] based on manual field measurements. In future investigations, it may be preferable to perform radiative transfer simulations using trees digitized directly from terrestrial laser scanning measurements [99,124,185,187] or digital photogrammetry [241] since these techniques can enable the importation of biologically realistic trees.

Precision of ALS Metrics
The DART simulations produced discrete-return height metrics that were consistent across 100 forest inventory plots, with little evidence of deleterious perturbations induced by varying acquisition parameters or ground topography. These findings supported the findings of previous investigations indicating that discrete-return heights are robust to moderate variation in survey conditions [69,72,219] and may therefore provide reliable predictor variables for the indirect retrieval of forest inventory measurements across different ALS surveys [70,242,243]. Conversely, discrete-return density, canopy cover and the proportion of multiple canopy returns were sensitive to variations in sensor altitude, laser scanning angle, pulse repetition frequency and transmitted pulse duration, indicating that these survey characteristics play an important role in governing the stability of discrete-return metrics obtained from different ALS surveys. These latter findings indicate that canopy cover and density metrics are only directly comparable for discrete-return data collected under similar survey conditions and are less reliable predictor variables for the retrieval of forest biometrics. Furthermore, the proportion of multiple canopy returns acquired using different survey parameters was highly sensitive to variations in laser pulse duration and scanning angle, indicating that ALS metrics derived from intermediate and last returns are not consistent between surveys and should be avoided where possible. Our DART simulations therefore corroborate similar findings reported in controlled ALS surveys by Morsdorf et al. [63], Naesset [66], Jakubowski et al. [61], Keränen et al. [67] and Sumnall et al. [244].

Impact of Acquisition Parameters
The DART simulations demonstrated that lower pulse repetition frequencies and wider scanning angles can cause a downward shift in discrete-return heights as a greater proportion of returns are recorded at crown edges rather than treetops. Discrete-return height percentiles declined by ∼0.5 m as laser scanning angles increased from 20 • to 30 • , indicating that wider scan angles (>20 • ) may be sub-optimal for forests where tree apices are conical with well-defined leaders. Our simulations also revealed that wider scan angles (>20 • ) may produce inaccurate canopy cover measurements due to a reduction in discrete-return density. These findings are supported by Holmgren et al. [106] who recorded a similar downward shift in the ALS height percentiles and perturbations to canopy cover measurements at wider scan angles. Holmgren et al. [106] found that these effects were most pronounced at the edges of the ALS swath. However, the magnitude of these effects was also dependent on stocking density and crown morphology, with dense stands of trees with elongated crowns experiencing a smaller downward shift in discrete-return height percentiles due to the occlusion of the lower forest canopy.
Our results contrasted with those of Chasmer et al. [74], Hopkinson [75] and Naesset [66] who reported an upward shift in canopy height distributions at lower pulse repetition frequencies.
The authors attributed this trend to an increase in peak pulse power at lower pulse frequencies, causing an upwards shift in discrete-return heights as pulses did not penetrate far into the forest canopy before triggering a return signal. We found little evidence of these pulse interactions in our DART simulations, potentially reflecting inaccuracies in tree morphology, model parameterization or experimental design. Our results indicate that any rise in discrete-return height distributions as a consequence of increased pulse energy was overshadowed by the failure to detect treetops due to a reduction in discrete-return density at lower pulse repetition frequencies.
Ascending platform altitudes were associated with an upwards shift in discrete-return heights, particularly for height percentiles in the upper canopy (≥ h 75 ). This upwards shift was primarily caused by the widening of laser beam diameters from 0.12-1.2 m as platform altitude increased from 400 to 1500 m (beam divergence = 0.3-0.8 mrad), resulting in a greater probability of pulses directly intercepting treetops as demonstrated empirically by Naesset [72] and Morsdorf et al. [63]. The widening of laser footprints at higher altitudes may also have been responsible for an increase in canopy cover and the proportion of multiple returns (Figure 8) as more biological material was intercepted in each laser footprint. These results contrasted with the empirical observations of Goodwin et al. [73] and Naesset [66] who recorded fewer multiple returns at higher platform altitudes due the attenuation of laser pulse energy across wider footprints. Despite this, our results highlighted the importance of selecting appropriate laser beam divergence thresholds when surveying forests at varying platform altitudes [245]. While the dispersal of laser beams at higher platform altitudes increases the probability of intersecting treetops with conical apices, pulses with narrower beam divergence (0.3-0.5 mrad) are often used for topographic mapping as narrow pulses penetrate closed-canopy forests to enable the acquisition of well-distributed ground returns [78]. Since ALS surveys are seldom undertaken specifically for forest inventory applications, further work is needed to investigate how beam divergence and platform altitude can be optimized to maximize the utility of ALS surveys for forest mensuration while simultaneously fulfilling the needs of other data users.
A potentially novel finding from our research is that longer pulse durations can trigger a decline in the proportion of multiple canopy returns as Gaussian waveform decomposition failed to resolve discrete-return signals from near-ground structures as half-pulse duration increased beyond 5 ns. It is unclear whether the observed decline in multiple canopy returns was caused by: (i) reduced slant range resolution at longer pulse durations, or (ii) by echo-detection errors within the Gaussian waveform decomposition algorithm. However, these findings imply that the ability to modulate transmitted pulse durations may be useful when surveying: (i) closed-canopy forests where the acquisition of reliable ground returns is problematic [211,233], or (i) multi-layered forests where understory vegetation is only partially detected [37,246]. Unfortunately, this is seldom feasible in practice since transmitted pulse durations are often sensor-specific and non-adjustable. More significantly, our findings indicate that Gaussian waveform decomposition may be sub-optimal for near-ground echo detection, lending support to the development of alternative peak-detection methods capable of extracting near-ground laser returns for more detailed measurement of understory vegetation and ground surface topography [247,248].

Conclusions
In this study, a variance-based global sensitivity analysis was undertaken to: (i) investigate how accurately the DART model was able to reproduce small-footprint ALS measurements acquired over commercial conifer plantations, and (ii) examine the impact of sensor and survey characteristics on the precision of discrete-return metrics commonly used for the derivation of forest inventory measurements. Our investigation demonstrated that DART could replicate small-footprint discrete-return metrics collected over 100 forest inventory plots dominated by Sitka spruce and Lodgepole pine with a high degree of accuracy and precision. These promising results demonstrated successful parameterization of the DART model for reproducing ALS data acquired by the RIEGL LMS-Q680i sensor in commercially important forests widespread throughout Britain and Ireland. Despite considerable modifications to survey settings and ground topography, the simulated discrete-return height distributions proved to be highly stable, indicating that discrete-return height metrics are potentially reliable predictor variables for the indirect retrieval of forest inventory biometrics across different ALS surveys. Conversely, ALS metrics such as discrete-return density, canopy cover and the proportion of multiple canopy returns were sensitive to variations in platform altitude, laser scanning angle, pulse repetition frequency and pulse duration. These latter findings indicate that canopy cover and density metrics are only comparable for discrete-return data collected under similar acquisition conditions, precluding their universal use for the derivation of forest inventory measurements across different ALS surveys.
Our study also highlighted some of the challenges and pitfalls of modeling complex forest canopies using synthetic three-dimensional trees generated from allometric relationships derived through manual field measurements. While our simulations replicated most ALS metrics with high fidelity, the proportion of multiple and ground returns was overestimated by 5.1% and 26.1% due to the misrepresentation of foliage. The accuracy of future simulations may therefore be improved by the importation of synthetic trees with biologically accurate shoot morphology, or the use of trees digitized directly from terrestrial laser scanning or photogrammetry. While our study successfully reproduced ALS measurements over structurally heterogeneous plantations, performing radiative transfer simulations over species-rich semi-natural forests remains a significant challenge for future research due to complexities associated with the replication of foliage morphology, density and orientation. Moreover, the parameterization of radiative transfer models such as DART requires considerable a priori knowledge of sensor specifications that is not always available for commercial ALS devices. Despite these caveats, we conclude that DART is a robust model for simulating discrete-return ALS measurements provided that realistic three-dimensional trees are imported into radiative transfer simulations. Funding: This research was funded through the Knowledge Economy Skills Scholarships (KESS2) initiative (kess2.ac.uk). KESS2 is a pan-Wales higher education initiative led by Bangor University and part-funded by the Welsh Government's European Social Fund (ESF) convergence programme for West Wales and the Valleys.

Acknowledgments:
We kindly thank the forest managers and inventory personnel at Coillte Teoranta for providing the ALS data and field measurements used to calibrate the DART model. We are also grateful to Jean-Philippe Gastellu-Etchegorry and Tiangang Yin for their advice and support in using the DART model. Finally, SuperComputing Wales (www.supercomputing.wales) and Colin Sauze are thanked for supporting the project through the provision of a High-Performance Computing (HPC) facility on which most of the DART simulations were performed.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

ALS
Airborne