The study was carried out on a 4.8-hectare cherry orchard located in the Provence-Alpes-Côte d'Azur region of southeastern France, centered at (44°11′34.8^′′N, 5°09′24.4^′′E), see Figure 1. The area has a Mediterranean climate, characterized by hot and dry summers with average daily temperatures around 30 °C, and cold winters with an average daily temperature below 10 °C. Annual precipitation ranges from 650 to 900 mm, with heavy rains usually occurring in September and October (Rouault et al., 2024). The orchard is organized into rows with a 7 m spacing and 5.5 m between trees, with wild grass covering the inter-row spaces. In 2022, the grass was mowed twice, late April and late May, and then left to dry out. Re-greening of the inter-row grass was observed starting in early September, followed by slow subsequent growth. The orchard contains five varieties of cherry: Prime-Giant, Belge, Belise, Folfer, and Summit, planted to provide early, mid, and late season yields to ensure a prolonged harvest period of approximately 40 days. The trees are planted at a density of 240 trees per hectare, with an average canopy dimension of 5.9 m (length) × 5.3 m (width) × 3.5 m (height). This corresponds to a tree density of 0.024 trees per square meter. The orchard is equipped with a drip irrigation system with a designed flow rate of 31.4 L per hour per tree. The irrigation system operated for 6 h per week in April, increased to between 1 and 3 h per day during May, June, July, and mid-August, and then was reduced to between 3 and 10 h per week, based on weather conditions and local authority regulations. A data logger was installed near the center of the orchard to collect in situ soil moisture data. Cherry trees were pruned and thinned in mid-November, and all work was completed by the end of November.

The CRNS used in this research is the CRS2000B from Hydroinnova LLC (2024). The system consists of two detector tubes, each approximately 1.5 m tall and filled with boron trifluoride gas for high neutron absorption. One detector is moderated by a 2.5 cm thick polyethylene layer (moderated detector), while the other remains unmoderated (bare detector). When neutrons enter the tubes, some are absorbed by the gas and generate electrical signals that are recorded. The bare detector measures low energy (thermal neutrons), while the moderated detector measures higher energy (epithermal neutrons). In February 2022, the CRNS was installed in the middle of the cherry orchard, mounted on a stand at 50 cm height. The bottoms of the detectors align with the bottom of the canopy, and the tops of the detectors are still below the average canopy height (around 3.5 m). A data logger records the total neutron count as well as with the average temperature, relative humidity, barometric pressure, and total rain data.

The raw neutron (N_raw) data are collected in hourly intervals and reported in units of counts per hour (cph). Before analysis, the raw data undergo a series of filtering and correction steps to remove invalid measurements and account for the influence of environmental variables. Outliers in N_raw data were removed following established CRNS data pre-processing protocols. Specifically, outliers (i.e., N_raw < 50 cph or N_raw>10, 000 cph) and values exceeding the 24-h rolling mean by more than twice the standard deviation were excluded (Bogena et al., 2022). After outlier removal, correction factors were applied following (Zreda et al., 2012):

where f_p, f_i, and f_h are the correction parameters for atmospheric pressure variations, incoming cosmic-ray neutron intensity variation (Zreda et al., 2012) and atmospheric water vapor (Rosolem et al., 2013), respectively. The correction applied for incoming cosmic-ray neutron intensity accounts for elevation and geomagnetic latitude using the scaling approach of McJannet and Desilets (2023).

Some studies suggest correcting thermal neutron counts only for atmospheric pressure and absolute humidity (Jakobi et al., 2018). Sensitivity analyses without the incoming neutron correction showed no major impact on the relationship between thermal neutron counts and vegetation metrics. These results are presented in the Supplementary Figures S8, S9. In addition, recent studies by Schrön et al. (2024) and Rasche et al. (2023) proposed that the water vapor correction factor for thermal neutrons is approximately 40% of the value used for epithermal neutrons. Therefore, following Schrön et al. (2024), the water vapor correction factor was adjusted from the standard correction factor (0.0054 m³g⁻¹) to the reduced factor (0.0021 m³g⁻¹) for thermal neutrons. The factors are calculated using the following equations:

where the coefficient α is:

where I (ctss⁻¹) is the neutron intensity measured by the Jungfraujoch (JUNG) monitoring station in Switzerland (available at www.nmdb.eu) at a given time, and I_ref (ctss⁻¹) is an arbitrary reference value used for normalizing cosmic-ray intensity data. In this study, I_ref was set to 159.4 (ctss⁻¹), the daily average measured at Jungfraujoch on 2011-05-01. h is the absolute air humidity (gcm⁻³) and h₀ (gcm⁻³) is the average air humidity over the study period. P is the measured air pressure (mbar) and P₀ is the average pressure over the study period. The neutron barometric coefficient, β, was estimated using a tool developed by CRNSLAB using the calculation described in McJannet and Desilets (2023) (available at crnslab.org/util/intensity.php) and has a value of 136 mbar^–1 at the study site. Lastly, daily mean N_epi and N_th values were calculated to compare to the UAV-LiDAR-derived PAI and satellite-derived vegetation indices. The daily means were computed from the hourly CRNS data, and the daily standard deviation is considered an indication of uncertainty.

Six HydraProbe-II sensors (Stevens Water Monitoring Systems Inc., Portland, Oregon, USA) were installed approximately 40 m from the CRNS and collected soil moisture data every 30 min. Their location is indicated by the orange rectangle in Figure 1. The sensors were divided into two sets, each consisting of three sensors installed at 5, 15, and 30 cm depths. One set was installed within the tree rows, to represent the irrigated portion of the field which accounts for 17% of the field area. The other set was installed in the inter-row spaces, which represents the remaining 83% (non-irrigated) field area.

The mean soil moisture value at each depth was calculated using area-based weights of 17 and 83% before calculating a daily mean. To account for depth-related contributions to the cosmic-ray neutron signal, the weighting approach proposed by Schrön et al. (2017) was applied using weights of 0.7, 0.25, and 0.05 for the 5, 15, and 30 cm layers, respectively.

Plant development in this study was primarily quantified using the plant area index (PAI), derived from UAV (Unmanned Aerial Vehicle)-based LiDAR (light detection and ranging) data. UAV-based LiDAR provides unique three dimensional canopy structure information that can be used to calculate PAI. PAI (m² m⁻²) is defined as the surface area of all above-ground plant components per unit of ground area. It includes both green and non-green leaves as well as woody parts of the plant, such as stems and branches (Jonckheere et al., 2004). PAI is widely used for vegetation monitoring and modeling purposes (Grau et al., 2017). UAV-LiDAR data are suitable for estimating tree structural properties such as height, and canopy volume and density (Farhan et al., 2024). However, data collection is labor-intensive and time-consuming. Although it does not provide a direct estimate of vegetation water content, variations in height, canopy volume and leaf density provide a measure of changes in above ground biomass and hence water content, which influences hydrogen and neutron counts. PAI is therefore used to characterize canopy development throughout the phenological cycle of the orchard.

UAV-based LiDAR data were acquired in 2022 with at least one acquisition per month. The UAV-LiDAR system used was a DJI Matrice 300 RTK quadcopter (DJI Technology Co., Ltd., 2020) equipped with a DJI L1 LiDAR sensor (DJI Technology Co., Ltd., 2021). The flights were carried out at an altitude of 50 m, with a velocity of 2.5 m s⁻¹ and a 50% side overlap to ensure complete coverage of the orchard and providing a point density of roughly 3,300 points m⁻². Each laser pulse emitted by the LiDAR sensor can generate up to three returns, reflecting from different layers of the canopy such as leaves, branches, or the ground. This multi-return capability helps to capture the vertical structure of the trees. Additionally, a real-time kinematic (RTK) system was utilized for georeferencing, supplemented with nine Ground Control Points (GCPs) to enhance positional accuracy. The UAV was also equipped with a DJI Zenmuse P1 RGB camera (DJI Technology Co., Ltd., 2019). The photos were processed with Agisoft Metashape software (Agisoft, 2024) for geometric rectification and photogrammetry, resulting in detailed orthomosaics of the orchard. The RGB images collected during the campaign capture the progression of tree phenological stages throughout the year, from leafless to full canopy development and back to leafless (see Figure 2). Figure 2 highlights the main transitions in canopy phenology, including dormancy (February–March), flowering (April), fruit growth (May–June), fruit ripening and harvest (July–August), and leaf senescence (October–December).

A voxel-based method was used to estimate PAI from the UAV-LiDAR data. The point cloud was processed using the open-source leafR package (version 0.3.5), which was developed to map canopy metrics, including PAI and height using a voxel-based approach (Almeida et al., 2021). The point cloud was represented in a grid of 3-dimentional cubes (voxels), and PAI was estimated for each voxel. PAI was estimated using an application of the Beer–Lambert law to relate canopy density to the attenuation of radiation through vegetation. For the discrete-return LiDAR data, cumulative PAI at a given canopy depth was derived by approximating light attenuation as the ratio of LiDAR pulses entering and exiting successive canopy layers (i.e., voxels), following Bouvier et al. (2015). To compute PAI in each voxel, the acquired point cloud data were structured into a 3D voxel grid (70 cm × 70 cm × 30 cm) using the x, y, z coordinates of each LiDAR return. The voxel size of 70 cm × 70 cm was chosen to maximize the number of LiDAR points used for PAI estimation given the average canopy dimension of approximately 5.3 m × 6 m. Then, for each column of voxels, the number of points within each voxel (from voxel 1 to voxel n) was calculated. PAI was estimated for each voxel (voxel i) as the ratio of the cumulative number of points from the top voxel down to voxel i (voxels 1 to i) to the cumulative number of points from voxel i + 1 to the bottom voxel (voxels i + 1 to n). The PAI for each column was then calculated as the sum of the individual voxel-based PAI values in that column. The outcome was a PAI map for the entire plot with a horizontal resolution of 70 cm × 70 cm. Previous studies have shown that 86% of the measured thermal neutrons were thermalized within a radial distance of 43–48 m from the CRNS (Jakobi et al., 2021). Therefore, voxel-based PAI values within a 45 m radius of the CRNS were averaged for comparison with thermal neutron counts.

Multispectral satellite data are also employed to provide additional insights into the orchard phenology and development. While satellite indices primarily capture signals from the top canopy surface, their broad spatial coverage and frequent revisit times make them highly suitable for vegetation monitoring. Satellite remote sensing data are widely used to derive vegetation indices such as normalized difference vegetation index (NDVI), normalized difference water index (NDWI), and normalized difference red edge (NDRE), which can be used to assess canopy condition and estimate biomass. The leaf area index (LAI), which is the ratio of the green leaf area to the ground area (in m² m⁻²), can be estimated from these indices.

A total of 164 cloud-free PlanetScope SuperDove surface reflectance images were collected during the UAV survey period between February 15 and December 1, 2022. These images were downloaded from the Planet Explorer platform (Planet Labs Inc., 2023) and clipped to a radius of 45 m around the CRNS to match the average footprint of the N_th detector (Jakobi et al., 2021, 2022). For days with multiple image acquisitions, the images were averaged, resulting in 98 images being retained for analysis. The data include eight spectral bands (i.e., red edge, red, green, green I, yellow, blue, coastal blue, and near-infrared) with a spatial resolution of 3 m.

In addition, the images were harmonized with Sentinel-2 using the “harmonize” imagery option available in Planet Explorer to ensure the spectral consistency of the vegetation indices (Kington and Collison, 2022; Houborg and McCabe, 2018). This approach corrects for the difference in surface reflectance between PlanetScope and Sentinel-2 common bands (Green, Red, Red Edge, and Near-Infrared), due to differences in sensor wavelength characteristics as provided in Table 1 (Baldin and Casella, 2024; Houborg and McCabe, 2016) .

Surface reflectance Level 2 Sentinel-2 imagery was obtained via the Google Earth Engine, following the correction procedure outlined in (Hagolle et al., 2008). Imagery was obtained from February to December 2022, with a temporal resolution ranging between 3 and 6 days. Due to cloud cover, a total of 76 suitable images were available during the 11-month study period. The resolution of Sentinel-2 images is 10 m for the green, red and NIR images, and 20 m for the red edge band (European Space Agency, 2020).

Four key vegetation indices (Normalized Difference Vegetation Index “NDVI,” Normalized Difference Red Edge Index “NDRE,” Normalized Difference Water Index “NDWI,” and Leaf Area Index “LAI”) were calculated from both PlanetScope and Sentinel-2 data. These indices (Equations 5–7) are computed using common spectral bands common to both sensors (Near-Infrared, Red, Green, and Red Edge), and are commonly used in vegetation monitoring. The NDVI is widely used to assess vegetation health and is calculated as follows:

where NIR is the near-infrared reflectance and Red is the red band reflectance. The NDWI relates to vegetation water content and is calculated as:

where Green is green band reflectance. NDWI has been shown to be more sensitive to canopy water content and biomass (Contreras et al., 2025). NDRE is sensitive to chlorophyll content and detects early plant stress before it is visible in NDVI, using a narrower red-edge waveband and calculated as:

where RedEdge is red-edge band reflectance. The use of these indices is based on their correlations with field-measured biomass in orchard systems (Panumonwatee et al., 2025).

Leaf Area Index (LAI) provides a measure of canopy density (Tian et al., 2023), and was estimated from Sentinel-2 and PlanetScope images using the neural network-based algorithm of biophysical variables (BVNet) (Weiss et al., 2002). This algorithm uses reflectances from three spectral bands B3 (Green), B4 (Red), and B8 (Near Infrared) and incorporates solar and sensor angles to refine LAI estimations (Courault et al., 2021). The algorithm utilizes a neural network that has been trained on simulated spectral reflectance, using the SAIL radiative transfer model as described in Weiss et al. (2002). Mukhtar et al. (2022) found that the LAI obtained through the BVNet algorithm was capable of tracking the dynamics of leaf development. Based on satisfactory evaluations of BVNet against ground-based measurements performed by Courault et al. (2021), BVNet has been incorporated into the ESA (European Space Agency) S2 toolbox (Abubakar et al., 2023).

Thermal neutron counts were correlated with each of the satellite data-based vegetation indices and the LiDAR-based PAI, and the goodness of fit was evaluated using the determination coefficient (R²). As demonstrated by Jakobi et al. (2021), approximately 45% of thermal neutrons originate from within a 5 m radius of the detector, with 86% of cumulative contributions coming from a radial footprint of about 45 m. Supplementary Figures S1–S4 show results obtained assuming smaller plot radii around the CRNS (10 and 25 m), and show that the correlation increases as the footprint is increased to this value. Therefore, results in the following section are based on a radial footprint of 45 m.

Representative maps of UAV LiDAR-derived PAI are shown in Figure 3 at both individual tree (panel i) and footprint scales (panel ii) to highlight the spatial variability in PAI across different times of the year. The images in Figure 3i visualizes the three-dimensional LiDAR point clouds for an individual tree at three key phenological stages. The X and Y-axes represent relative horizontal coordinates in meters, used to illustrate the spatial extent and dimensions of the tree rather than absolute geographic position. The Z-axis represents canopy height in meters. On 2022-03-29 (dormancy), the tree consists of a sparse structure with minimal foliage. By 2022-05-30 (ripening period), the canopy has developed into a dense structure with nearly maximum leaf coverage. On 2022-12-01 (post-senescence), the tree returns to a bare structure with only woody elements. PAI maps at the CRNS footprint scale are shown in Figure 3ii with X and Y-axes representing the relative coordinates in meters (X = east-west, Y = north-south), with the CRNS located at the center. On 2022-05-30 (Figure 3ii-a) when leaves are fully developed and fruit is approximately halfway through ripening, PAI values are high across the orchard. Individual trees are discernible, with peak values up to 5.8 m² m⁻² near the canopy center and decreasing PAI with distance from the trunk due to the decrease in branch density. In contrast, on 2022-12-01 (Figure 3ii-b) PAI is uniformly low < 0.7 m²m⁻² after senescence and pruning. At this stage, the PAI captures the woody structural elements of the trees. In both maps, the white areas indicate zones of “no data” where the PAI of the canopy is zero due to gaps in the canopy cover. In these zones, the short inter-row vegetation (grass) is located entirely within the bottom voxel making it impossible to obtain a reliable estimate of PAI.

Figure 4 shows the in-situ rainfall and soil moisture, observed PAI and neutron counts during the study period. Soil moisture (Figure 4a) exhibits some seasonal variation, with values around 0.30 m³ m⁻³ during the winter and spring months (January–April), followed by a gradual decline to approximately 0.10–0.15 m³ m⁻³ during the summer period. Several rainfall events occurred throughout the year (indicated by gray bars). The larger rainfall events produce soil moisture drydowns in April, July, and September. Vegetation growth stage is indicated in Figure 4a. The PAI in the 45 m CRNS footprint is lowest during dormancy (0.5 m² m⁻²) and increases from early spring (March) to reach its peak during early summer (end of June) (2.8 m² m⁻²) when the canopy is fully grown ready for harvest. After harvest, between July 1st and mid-August, the average PAI of the footprint decreases in response to the 50% reduction in irrigation. Irrigation was fully stopped in mid-September, resulting in leaf rolling and a general canopy shrinkage and a decline in PAI values (Rouault et al., 2024). Between mid-August and end-September, PAI values increased as a result of multiple rain events that eased the water stress, rehydrating the canopy and allowing leaf expansion. The onset of senescence, pruning activities, and natural leaf loss during November and December resulted in a sharp decline in PAI from early October.

Figure 4b shows the difference in response for thermal and epithermal neutrons to the changes in soil moisture and vegetation development. The epithermal neutron counts exhibit a consistent and proportional response to soil moisture dynamics, with sharp dips in response to precipitation followed by a slow recovery during drydowns. In contrast, sharp dips in N_th are only apparent from September onwards. During three rainfall periods between April and July, thermal neutron counts do not show a clear decrease following the increase in soil moisture due to the rain. N_th follows the PAI more closely, increasing from dormancy through flowering and fruit set to peak canopy development. Between July and September, as the soil dries, epithermal neutron counts increase in response to decreasing soil moisture, whereas thermal neutron counts and PAI are both decreasing.

In November, thermal neutron, like epithermal neutron counts, decrease following heavy rainfall as a response to a strong increase in soil moisture. The drops in N_th values occur due to the increased presence of hydrogen in the recently wetted soil and intercepted water in the canopy (i.e., water droplets on leaf surfaces). While the amount of water intercepted by vegetation is typically small (i.e., less than 2 mm) and short-lived, water puddles formed after rainfall due to soil saturation have a more significant impact by enhancing thermal neutron absorption and reducing measured count rates (Schrön et al., 2017). This leads to a temporary decorrelation between the PAI and thermal neutron measurements. After that, while epithermal neutron counts stabilize under relatively constant soil moisture conditions, thermal neutron counts continue to decrease in response to leaf falling season. Overall, during the main growing season, thermal neutron variability seems to be primarily driven by vegetation dynamics rather than soil moisture, closely following changes in PAI and visually observed indications of water stress (leaf rolling and expansion). Figure 4 shows that epithermal neutron counts respond to soil moisture changes throughout the season. Thermal neutrons only exhibit a clear response to soil moisture dynamics following significant rainfall events and during the senescence stage.

The relation of thermal and epithermal neutrons to soil moisture and vegetation are further examined in Figure 5, which presents the correlation between daily neutron counts and in-situ SM measurements. Figure 5a shows that N_epi is strongly correlated with in situ SM (r = −0.95, R² = 0.91), with data points mostly clustered along the regression line. The high R² value indicates that approximately 91% of the variance in N_epi can be explained by soil moisture alone, with limited influence of other hydrogen sources such water within the vegetation. This strong relationship explains why N_epi was very responsive to SM variations throughout the year, as expected from CRNS theory and existing literature.

Figure 5b shows that N_th has a weaker relationship with in situ SM (r = −0.78, R² = 0.61). SM clearly influences N_th. However, as shown in Figure 4, the response to soil moisture occurred primarily in response to very large rainfall events and during senescence. This explains why R² between N_th and SM is lower than R² between N_epi and SM.

The data points are scattered around the regression line, suggesting that other hydrogen pools are contributing to the variance. For example, in the soil moisture range (0.23–0.27 m³ m⁻³) N_th varies from 1,070 to 1,140 cph for similar soil moisture levels. At SM values between 0.1 and 0.2 m³ m⁻³, N_th varies from approximately 1,150 to 1,240 cph, a range of 90 cph. These observations suggest that variations in N_th during the growing season are primarily driven by something other than soil moisture fluctuations.

Figure 6 shows the relationship between PAI and neutron counts, using color to indicate the time of year. Figure 6a shows that N_th exhibit a strong linear correlation with PAI (R² = 0.86). Increasing PAI indicates an increase in leaf and/or stem area (and volume), thereby an increase in hydrogen stored as water in the vegetation. The slope indicates that N_th is sensitive to this change. Figure 6b shows that the relationship between N_epi and PAI is weaker (R² = 0.29). This is consistent with existing literature, and the previous figure which showed that N_epi was primarily sensitive to soil moisture.

Note that the PAI data were resampled to 3 m resolution to match PlanetScope imagery. However, at 3 m resolution, each voxel includes both vegetation and empty ground space, which leads to a considerable number of ground returns and hence lowers the PAI estimation. Supplementary Figures S5–S7 show that the relation with thermal neutron counts remain almost unchanged regardless of resolution, and the correlation between PAI and N_th at 3 m voxel resolution is still high (R² = 0.84).

Following the example of Brogi et al. (2025), a LOOCV analysis was performed to test how well PAI could be predicted using thermal neutron count rates. The cross-validation yielded an R² of 0.73 and an RMSE of 0.4 m² m⁻². This RMSE corresponds to about 14% of the observed range in PAI values (see Figure 4). The strong correlation between N_th and PAI suggests that N_th may be suitable for monitoring plant development. However, it remains important to acknowledge the potentially confounding influence of soil moisture, as reported in Figure 4.

The vegetation index maps showed patterns consistent with the LiDAR-PAI maps, albeit at coarser resolution. As an example, Figure 7 shows NDRE maps from PlanetScope and Sentinel-2 at two distinct periods. PlanetScope-derived NDRE in Figure 7a shows a spatial distribution similar to the PAI map shown in Figure 3a with two distinct clusters of values, with lower values in the lower right quadrant. This is likely due to missing trees or recent replanting (smaller canopy size), which results in more inter-row (bare or vegetated soil) being captured by the satellite. Depending on the inter-row conditions, this surface can have lower or higher NDRE values compared to the trees. For instance, in Figure 7c, the pixels in the lower right quadrant show higher NDRE values (approximately 0.45) compared to the rest of the plot (which is around 0.3). This is likely due to the presence of grass in that quadrant, which contributes more to the vegetation indices by reflecting more in the near-infrared band compared to trees in the senescence stage (see Figure 2). Similar patterns are observed in Sentinel-2 images shown in Figures 7b, d despite their lower resolution.

Figure 8 compares the temporal variation of N_th and vegetation indices (NDVI, NDWI, NDRE, and LAI) averaged over the 45 m radius area of the CRNS footprint. The PlanetScope and Sentinel-2 vegetation indices capture the seasonal dynamics in the orchard, such as the stages of dormancy, growth and senescence. The dip between June and mid-August was caused by leaf rolling and stress, followed by a recovery between August and October after multiple rainfall events and a steady decline starting in October with the onset of leaf loss (Rouault et al., 2024).

For example, the NDVI shows a marked increase from 0.5 in February to 0.8 by June which reflects vegetation development during this period. This is followed by a decline to 0.63 between June and August, likely due to water stress and leaf rolling. Between August and October, the NDVI increases again, reaching values between 0.64 and 0.8, driven by rainfall and recovery. However, it decreases steadily afterward, aligning with the onset of leaf loss. PlanetScope and Sentinel-2 show globally similar trends for the various calculated vegetation indices, despite small differences in absolute values, most noticeably for LAI.

The temporal variations of the vegetation indices reveal sensitivity to environmental conditions. For example, on 26 September in Figure 8, a sudden drop in vegetation indices was observed after heavy rain in the field. Such a drop may be explained by the temporary presence of water droplets on the leaf surfaces after rainfall. When water covers the leaves, it forms a thin layer that enhances the absorption in the near-infrared band before it can be reflected, while slightly affecting the red and green bands. That reflection mechanism causes a decrease in vegetation indices based on high near-infrared and low visible reflectance (Gao et al., 2025). Differences between the indices reflect their differing sensitivity to, e.g., chlorophyll, water content, etc. Although similar dynamics are observed for the PlanetScope and Sentinel-2 images, PlanetScope provides finer spatial details and greater temporal density. As with the LiDAR-derived PAI, Figure 8 shows that N_th covaries with vegetation indices, with a decorrelation in the trends observed when N_th responds to rain events.

To allow a direct comparison, the correlation between multispectral-derived vegetation indices and N_th was calculated using data from the same dates when the PAI was collected, see colored points in Figure 9. Weaker correlations with N_th, compared to PAI, were obtained with LAI derived from Sentinel-2 (R² = 0.50) and PlanetScope (R² = 0.59), see Figure 9a. This was expected because optical reflectance primarily captures the top surface layer of the canopy, whereas the entire canopy influences the neutron count. The correlation between N_th and vegetation indices is likely because the upper canopy is representative of the whole canopy. The purpose of these comparisons is to evaluate the relationship between N_th commonly used vegetation indices. While N_th has shown strong correlation with PAI in this orchard system, vegetation indices from satellite data have advantages in terms of spatial coverage and temporal frequency. Therefore, understanding how well vegetation indices correlate with N_th helps assess whether N_th could potentially complement or substitute for vegetation indices in monitoring orchard development, particularly when satellite data are limited by cloud cover or acquisition frequency. Among all of the vegetation indices examined herein, PlanetScope-derived NDRE, see Figure 9d and NDVI, Figure 9b showed the highest correlation with N_th, with R² = 0.69 and R² = 0.61, respectively. However, a lower correlation was observed between N_th and NDWI with R² = 0.48, see Figure 9c. PlanetScope-derived vegetation indices show consistently higher correlations with N_th compared to Sentinel-2. This is most likely due to PlanetScope's finer spatial resolution, which captures more pure pixels of vegetation compared to Sentinel-2's 10-m resolution that is more affected by mixed soil–canopy pixels.

The stronger correlation with NDRE compared to other indices is attributed to the sensitivity of the red-edge band, used in the NDRE calculation, to chlorophyll content and canopy structural characteristics. A recent study reported that the red-edge bands show the highest sensitivity (>85%) for structural leaf-parameters among multiple spectral bands (Yan et al., 2025). In contrast, Sentinel-2 LAI, NDVI, NDWI, and NDRE show similar, lower correlations with N_th, likely due to the uncertainties of model-based retrievals in heterogeneous orchard environments.

While the correlations presented in Figure 9 are based on dates coinciding with UAV-LiDAR acquisitions (colored points), the full temporal dataset (gray points) provides additional information about the N_th-vegetation index relationships. For example, the scatter in the complete dataset shows potential nonlinear, possibly saturation effects at higher vegetation index values, clearly observed in the NDVI subplot in Figure 9b where NDVI values plateau around 0.7–0.85 even as N_th continues to increase beyond 1,150 cph. Similar effects are observed in NDWI, as shown in Figure 9c (plateauing around 0.65–0.80), and NDRE (around 0.50–0.60). Optical indices become less sensitive to canopy changes at high biomass levels. Additionally, some of the temporal decorrelation between N_th and vegetation indices occurs during and immediately after rainfall events (as discussed in relation to Figure 8), which contributes to the scatter in the full dataset. The lower correlation between N_th and the multispectral-derived vegetation indices compared to the LiDAR-derived PAI is due to the fact that N_th is affected by the entire canopy while optical reflectances are mainly influenced by the top surface of the canopy.

Results in Figures 4, 5 indicate that while there is a high correlation between N_th and PAI, there is also sensitivity to soil moisture. In order to use N_th directly to predict PAI, a correction would be needed to account for the effect of soil moisture on N_th. One solution would be to establish a correction term based on a simple fit between N_th and SM during a period of low PAI. Such a correction term could be applied to the time series to obtain a SM-corrected thermal signal with which to predict PAI.

Obtaining such a correction from the observations in this study is challenging due to covariance between PAI and SM (SM tends to be high when PAI is low and vice versa) and the limited SM variability during the low PAI period (See Figure 4). During this study, SM variations occurred in the period between mid-April and October, the period in which PAI and the vegetation indices are also most dynamic (Figures 4, 8). Extending the study duration provides no guarantee that soil moisture variability would be observed during a low PAI period. Therefore, relying on calibration of a simple, empirical correction term is not always practical.

Figures 5, 6 show that both N_th and N_epi depend on both SM and vegetation, which is consistent with previous studies focused on soil moisture and irrigation estimation (Jakobi et al., 2022; Brogi et al., 2023; Li et al., 2019; Baatz et al., 2015; Brogi et al., 2025; Köhli, 2026). N_epi is more strongly correlated with SM than with vegetation (PAI and vegetation indices in this study). Hence, N_epi is used for SM estimation, acknowledging that the effect of vegetation must be taken into account (Morris et al., 2024). This study showed that N_th is more strongly correlated with PAI and vegetation indices than with SM. Therefore, it could be useful for vegetation monitoring provided the influence of soil moisture can be taken into account.

The use of a simple, empirical correction term may be challenging due to the seasonal variation in soil moisture and its covariance with PAI. More sophisticated approaches, e.g., based on data assimilation, should be considered to exploit the complementarity of the N_th and N_epi observations.