Using repeat airborne LiDAR to map the growth of individual oil palms in Malaysian Borneo during the 2015 – 16 El Ni ˜ no

Making oil palm agriculture as efficient as possible is essential to ensuring that this economically important crop can be grown sustainably. To determine how oil palm growth rates vary across tropical landscapes, we used repeat airborne LiDAR data to map the height growth of > 500,000 oil palms in Malaysian Borneo over a two- year period coinciding with the 2015 – 16 El Ni ˜ no drought. Despite uncharacteristically dry and hot weather conditions, oils palms grew an average of 1.6 m yr (cid:0) 1 in height over this period. However, oil palm growth rates varied across the landscape and in relation to plant age, tending to be fastest for younger individuals growing closer to forest edges, further from rivers and at higher elevations. Our results highlight the ability of oil palms to grow even during periods of drought and showcase how cutting-edge remote sensing technologies can help improve the efficiency and sustainability of oil palm agriculture.


Introduction
African oil palm has earned some notable accoladesincluding 'tree of life', 'wonder crop' and 'natures gift to man' -as almost every part of the plant has a use (Gillespie, 2012;Okolo et al., 2019). Oil palms produce ten times more oil per hectare than their closest competitor, soybean (Mayes, 2020;Sumathi et al., 2008), making oil palm a key crop worldwide for the production of vegetable oil over the last century (Okolo et al., 2019). Nowhere is this truer than in parts of SE Asia (Brad et al., 2015), for instance Indonesia-where palm oil accounts for around 10% of total national exports and is central to the livelihoods of countless people (Beyer et al., 2020). However, this rapid rise in the demand for oil palm has come at a huge environmental cost, with the growth of the oil palm industry leading to the large-scale conversion of regions of tropical forest into monoculture plantations. On the island of Borneo, over 3 million ha of primary forest were felled to clear space for oil palm agriculture between 2000 and 2017, with substantial impacts on biodiversity and carbon storage (Beyer et al., 2020;Dislich et al., 2017;Gaveau et al., 2016;Koh and Wilcove, 2007;Qaim et al., 2020). It has been evidenced that oil palm plantations hold less than 50% of the vertebrate species of primary forests and hold far lower overall biodiversity than both primary and degraded secondary forests (Fitzherbert et al., 2008). This puts many of the ecosystem services provided by stable and diverse tropical landscapes at risk, compromising soil retention and nutrient cycling with adverse long-term impacts on the productivity of the oil palm plantations themselves (Dhandapani et al., 2019;Dislich et al., 2017). Moreover, intact tropical rainforests can hold as much as 270 tonnes more carbon per hectare than oil palm plantations (Katayama et al., 2013;Khoon Kho et al., 2015;Sayer et al., 2012). Consequently, the continued expansion of oil palm agriculture is weakening the tropical forest carbon sink, which has performed a crucial role in slowing the pace of climate change over the past half-century (Qaim et al., 2020).
The solution to this challenge is not to simply ban or replace palm oil. Not only would this have devastating socio-economic impacts, but the environmental impacts of vegetable oil production would simply be displaced elsewhere. Given how productive oil palm is compared to alternative crops, this could end up doing more harm than good for both climate and biodiversity (Beyer et al., 2020). Instead, it is vital to find ways to produce oil palm more sustainably. One way to do this is to increase efficiency by maximizing yields within existing oil palm landscapesthe so-called 'land sparing' approach (Grass et al., 2019;Phalan et al., 2011). Research shows that there are often significant gaps between actual and potential oil palm yields. Harvested yields rarely exceed 3 t oil ha − 1 yr − 1 , when estimates suggest that better cultivation practices could almost triple average yields to around 8 t oil ha − 1 yr − 1 (Corley and Tinker, 2015;Qaim et al., 2020;Woittiez et al., 2017). Closing these yield gaps would be a success for not just people, but the environment as well, since it would allow more oil to be produced from less land. However, while some of the large-scale drivers of these yield gaps are known, such as climate and soil fertility (Ajeng et al., 2020;Sarkar et al., 2020), far less has been established on what causes yields to vary locally within landscapesthe scale at which farmers can actually intervene by adapting planting strategies. For example, local factors such as topography can have a profound influence on soil water and nutrient availability, as well as on microclimatic components, including solar radiation, air temperature, and humidity-all of which directly constrain plant growth (Jucker et al., 2018c). Moreover, these same landscape features can exacerbate or dampen the effects of extreme climate events, such as the unusually hot and dry conditions associated with El Niño (Nunes et al., 2021).
A main challenge of identifying oil palm yield gaps across tropical landscapes is the scale of the plantations, which is often vast. This makes tracking the growth rates of oil palms from the ground both logistically challenging and prohibitively expensive (Mayes, 2020). One solution to this challenge is to leverage remote sensing technologies to map oil palm plantations from above, an approach that is becoming increasingly popular under the banner of precision agriculture (Mulla, 2013;Sishodia et al., 2020). One particularly promising technology in this regard is airborne LiDAR, which can capture the 3D structure of vegetation and the underlying terrain in great detail (Mulla, 2013). This makes airborne LiDAR an ideal tool for measuring vegetation height and biomass at large scales and exploring how features constrain their variation (Asner and Mascaro, 2014;Jucker et al., 2018bJucker et al., , 2018a. LiDAR data has also been used to automatically identify and measure the size of the crowns of individual plants, including oil palms (Dalponte and Coomes, 2016;Nunes et al., 2017). All of this suggests that by using repeat LiDAR surveys conducted at two or more points in time, it should be possible to track the growth of individual oil palms over time and explore how and why it varies across entire landscapes (Nunes et al., 2021).
To test this idea, we acquired repeat LiDAR data from a region in Malaysian Borneo, predominantly covered by oil palm, and mapped the height growth of >500,000 individual oil palms between 2014 and 2016. This period coincided with the global 2015-16 El Niño event (Jiménez-Muñoz et al., 2016), which resulted in unusually hot and dry conditions lasting multiple months across the region (Doughty et al., 2021;Nunes et al., 2021). Using these data, we addressed two main questions about oil palm growth and its variation across tropical landscapes: (1) How much did oil palms grow during this period characterized by unusually hot and dry conditions, and how variable were growth rates across the landscape? And (2) what developmental, ecological and landscape features drive variation in the growth rate of oil palms during this two-year period?

Study area
The data were acquired in conjunction with the Stability of Altered Forest Ecosystems (SAFE) project, situated in the Malaysian state of Sabah, in Borneo (Ewers et al., 2011). The SAFE project is one of the biggest ecological experiments in the world on ecosystem and biodiversity changes across tropical forests as a result of human modification, forest degradation, and fragmentation. The region has a tropical climate, with mean annual temperatures around 26 • C and annual rainfall between 2,600-3,000 mm (Jucker et al., 2018a). The SAFE landscape is highly fragmented and comprised of a variety of land-use types, including oil palm plantations of various ages, logged and fragmented secondary forests, and unlogged old-growth forests. Our study focuses on oil palm plantations in this region spanning an area of approximately 94 km 2 .
The site was affected by the 2015-16 El Niño drought, one of the strongest on record in Borneo (Doughty et al., 2021;Nunes et al., 2021) and many other tropical regions including Amazonia (Jiménez-Muñoz et al., 2016). In particular, the 2015-16 El Nino was anomalous not just because of lower-than-average rainfall, but also because of the sustained period of high temperatures (Doughty et al., 2021;Nunes et al., 2021). Combined together, these hot and dry conditions put forests in the region under intense and sustained physiological stress, leading to increased rates of tree mortality and leaf shedding (Doughty et al., 2021;Jiménez-Muñoz et al., 2016;Nunes et al., 2021).

2014 LiDAR data
LiDAR data were first acquired across the SAFE landscape in November of 2014 with a Leica ALS50-II LiDAR sensor flown by NERC's Airborne Research Facility. Data were collected as a discretized point cloud, with up to 4 returns per laser pulse (Jucker et al., 2018c). LiDAR point density varied considerably across the landscape, being highest over forested areas and lowest over oil palm plantations, which had a mean of 2.8 points m − 2 . For the purposes of our study, we classified point cloud data into ground and non-ground returns with the software 'LAStools' (https://rapidlasso.com/lastools). We then fit a digital elevation model (DEM) to the ground returns to generate a 1 m resolution raster before subtracting the DEM values from the non-ground returns to generate a normalized point cloud. Next, we generated a 0.5 m resolution canopy height model (CHM) from the normalised point cloud using the pit-free algorithm described in (Khosravipour et al., 2014). Further details on the acquisition and processing of these data can be found in (Jucker et al., 2018a). The DEM, CHM and point density data are archived online at: 10.5281/zenodo.4020697. In addition to the LiDAR data, true-colour RGB imagery was also acquired across the SAFE landscape using a Phase One iXU-RS 1000 100 MP digital camera mounted alongside the LiDAR scanner. Individual images were subsequently georeferenced, orthorectified, and stitched together into a mosaic spanning the same area as the CHM and DEM.

2016 LiDAR data
Our second LiDAR dataset was collected by the Global Airborne Observatory (GAO; formerly the Carnegie Airborne Observatory; Asner et al., 2012) in April 2016 as part of a larger-scale project which mapped aboveground carbon stocks across the state of Sabah (Asner et al., 2018). These flights were conducted at a higher altitude than in 2014, as the goal of the project was to maximize coverage of the region to best capture spatial variation in forest structure and aboveground biomass. However, the LiDAR point density over the oil palm areas considered in this study was almost identical to that of the 2014 data (3.3 points m − 2 on average), allowing robust estimation of canopy height changes across the study area. The processing of the point cloud data followed a similar approach to the 2014 data described above. We used LAStools to classify the point clouds into ground and non-ground returns, following which a 2 m resolution DEM and CHM were created. Further details of the processing and acquisition of these data can be found in (Asner et al., 2018).

Data processing
All subsequent data processing and analysis were carried out using a combination of QGIS (QGIS Development Team, 2016) and R (R Core Team, 2020). For a general overview of the workflow described below, see Fig. 1. First, to minimize errors due to misalignment between the two datasets, we used the Georeferencer plug-in in QGIS to manually align the 2016 CHM to the 2014 data. Following this, we cropped the full extent of the 2014 and 2016 CHMs so that only overlapping areas covering oil palm plantations were retained for further analysis. This was achieved using a shapefile of the SAFE project landscape marking the boundary of oil palm plantations, as well as creating a shapefile marking the area of overlap between the two datasets, as the 2016 data only covered a portion of the area flown in 2014 (Fig. 2). Additionally, areas of riparian forest surrounding rivers and areas affected by cloud cover were manually delineated and removed from the CHMs.

Crown segmentation algorithm validation
To develop an effective approach for identifying and delineating individual oil palms from LiDAR across the entire SAFE landscape, we created a validation dataset that would allow a comparison of the accuracy of different crown segmentation routines. An area of approximately 4.5 ha straddling the boundary between a mature and young oil palm plantation was chosen and by using both the RGB imagery and 2014 CHM data, we manually delineated all oil palm crowns within this area. A 5 m buffer was then applied to the area and we removed oil palms that fell within this buffer to avoid including individuals with crowns falling partially outside the extent of the imagery. This left a total of 410 manually delineated oil palms to assess the accuracy of the crown delineation algorithm (Fig. 3).
To automatically identify and delineate the crowns of individual oil palms in the 2014 CHM we used inbuilt routines in the lidR package in R (Roussel et al., 2020). This involved first using a local maximum filter (LMF) algorithm to locate the tops of individual palms and then applying a segmentation algorithm to delineate the border of their crowns in the CHM. For the purposes of crown segmentation, we compared five alternative algorithms that have been proposed in the literature: dal-ponte2016, li2012, silva2016 and watershed (as implemented in the lidR package) and the original version of dalponte2016 implemented in the itcSegments package. Based on this, we chose the dalponte2016 crown segmentation algorithm from the lidR package for all subsequent analyses, as it ranked first in terms of both accuracy and speed (see Table S1 in Supporting Information for details). This algorithm has been used successfully in the past to map individual oil palms using LiDAR (Nunes et al., 2017), and has been shown to generally outperform other segmentation routines applied to rasterized CHM data (Aubry-Kientz et al., 2019;Eysn et al., 2015). It also has the key advantage over approaches that work on point cloud data that the algorithm is computationally scalable across large areas such as those considered in this study (~100 km 2 ).
The segmentation routine is described in detail in Dalponte and Coomes (2016). Briefly, the LMF algorithm allows the user to specify a window size across which to search for local maxima (i.e., palm crowns), which we varied between 6 and 10 m. We used a fixed-size window with this range of values (as opposed to a variable-size window), as on average oil palms were observed to be spaced approximately 7-9 m apart regardless of their size, reflecting the typical planting configuration of oil palm plantations. Once the local maxima have been identified,  Crowns falling within a 5 m buffer from the edge of this area were excluded (purple polygons in c), leaving a total of 410 manually delineated crowns for validating algorithms (yellow polygons in c). The accuracy of the automatic oil palm segmentation routine using a window size of 9 m is shown in (d). Correctly segmented oil palms represent those where a single individual was found within a manually delineated polygon (yellow polygons in d), over segmented palms are ones where more than one individual was found within a manually delineated polygon (light green polygons in d), and omitted palms are ones where no individual was found within a manually delineated polygon (blue polygons in d). For additional details on the segmentation routine see Table S1 and S2 in Supporting Information. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the segmentation algorithm then uses a decision tree approach to grow a crown around each individual treetop (Dalponte and Coomes, 2016).
We assessed the accuracy of the segmentation routines by comparing the computer segmented crowns to those manually delineated by hand and calculating (i) the detection rate (i.e., the number of correctly segmented palms), (ii) the number of omitted palms (i.e., those which the algorithm failed to detect) and (iii) the number of over segmented palms (i.e., those which the algorithm incorrectly split into two or more palms) (Fig. 3). The primary goal was to keep the number of oversegmented palms as low as possible, as this would otherwise introduce a source of pseudoreplication in subsequent analyses. We therefore set an upper threshold of 2% for over segmented palms. This threshold of 2% was chosen arbitrarily as a compromise between minimising oversegmentation while keeping omission rates as low as possible.

Calculating height growth of individual oil palms across the SAFE landscape
The best performing segmentation routine was used to automatically identify individual oil palm crowns across the entire 2014 CHM. To calculate height change between 2014 and 2016, we resampled the 2014 CHM to 2 m to match the resolution of the 2016 data. The polygons of the individual oil palm crowns were then overlayed onto the CHMs from 2014 and 2016 and we extracted the maximum height value of pixels falling within each polygon. Finally, to calculate the height change (in m yr − 1 ) of each oil palm we subtracted the height in 2014 from that of 2016 and divided this by the time interval between surveys (17 months).
As part of this process, we excluded any palms that did not have any data attributed it for any of the features we investigated, before filtering palms that were <2 m , and with a crown area <4 m 2 in 2014 from the analysis. To filter out remnant palms left within the landscapes that could have been mistakenly counted as oil palms, any individual >18 m in height was also excluded (Barcelos et al., 2015). Finally, we used Mahalanobis distance as implemented in the Out-lierDetection package in R to identify and remove palms that exhibited an unrealistic change in height between 2014 and 2016 (see Fig. S1). These outliers could have been the result of a mismatch in the data, represent non-palm vegetation, or reflect palms that died or were severely damaged between surveys. In total, 14,144 palms were removed as outliers (2.1% of the total). We note that previous work has shown that LiDAR-based estimates of oil palm height are closely corrected to fieldbased measurements (Nunes et al., 2017), indicating that it should be possible to robustly estimates rates of oil palm height growth using repeat LiDAR acquisitions.

Drivers of local-scale variation in oil palm height growth
To explore the factors that could explain variation in oil palm growth rates across the landscape, we assembled data on a range of features linked to topography, ecological context, and plant development stage. First, the 2014 DEM was used to calculate several terrain metrics that have been shown to capture variation in hydrology, soil water availability, nutrients, local microclimate, and exposure to sun and wind (Jucker et al., 2018c). These included terrain elevation, slope, terrain ruggedness index (TRI), terrain position index (TPI) and aspect (see Figs S5-S10 in Supporting Information for maps of each of these spatial layers). Prior to calculating these metrics, we aggregated the DEM to a resolution of 10 m to smooth out any local artifacts and speed up subsequent computations. These predictors were derived using the raster package in R. Following this step, aspect values were cosine transformed to obtain a variable ranging in value between − 1 (corresponding to north-facing slopes) and 1 (south-facing) (Blonder et al., 2018).
In addition to these topographic metrics, we used a shapefile of rivers across SAFE to calculate the distance from rivers for all oil palms using the sf package in R (Pebesma, 2018). Similarly, the distance to the closest forest edge from each oil palm was calculated to test whether changes in microclimate, soil structure and pathogen loads related to proximity to intact forests might affect oil palm growth. Finally, to further test the effects of planting density and configuration, we determined gap fraction by aggregating the 2014 CHM to 20 m resolution and defining the proportion of pixels <1 m in height within a 20 m area around palms. These values were then subtracted from 1 to create a measure of 'canopy cover', where high values of canopy cover indicate palms surrounded by a high density of neighbours.
Multiple regression was used to determine how well each of these predictors contributed to explaining variation in oil palm growth across the landscape. To avoid issues with multicollinearity, we calculated Pearson's correlation coefficients (ρ) between all model predictors prior to model fitting (Fig. S2). Any that exceeded ρ > ± 0.6 were excluded from the analysis. Based on this, the following predictors were retained in the regression model: Palm height in 2014 (hereon in referred to as 'initial height'), canopy cover, distance to the forest edge, distance to rivers, ground elevation, TPI, slope, and aspect. All of the model predictors were scaled to have a mean of 0 and a standard deviation of 1 to allow their effect sizes to be directly comparable. As a post-hoc test of collinearity among model predictors, the variance inflation factors for the fitted model were calculated and confirmed as <2 for all predictor variables. We also tested whether individuals found in closer proximity to each other were more likely to exhibit similar growth rates than what may be assumed by chance, but found no evidence of spatial autocorrelation in the model residuals (Fig. S3).

Accuracy of oil palm segmentation
For the individual oil palm segmentation, a window size of 9 m proved the best compromise between minimizing over-segmentation while also ensuring as few individuals as possible were omitted by the algorithm (Fig. 3 and Table S2). With this configuration, the dal-ponte2016 segmentation algorithm implemented in the lidR package correctly segmented 67.1% of manually delineated oil palms and had an over-segmentation rate of just 1.7%, both of which were better than all other alternative segmentation algorithms we tested (Table S1). The segmentation accuracy of the dalponte2016 algorithm was greater for small palms <6 m in height (87.1%), whereas for mature individuals it was lower, at 57.1%. Reducing the window size to 8 m resulted in slightly higher segmentation accuracy (70.7%), but this came at the cost of doubling the rate of over-segmentation to >4% (Table S2).

Oil palm height growth and its variation
After applying the data quality filters described in the Methods, we delineated a total of 543,145 oil palms across the SAFE landscape and measured their height growth (Fig. 4). There was a significant difference in the mean initial height of palms (7.1 m) and their height in 2016 (9.2 m) (t = -254.08, P < 0.0001), with palms growing an average of 1.6 m yr − 1 between the two LiDAR flights (Fig. 5). However, there was considerable variability in the rate of height growth across the study area, with 90% of values ranging between 0.5 m yr − 1 (5th percentile) to 2.7 m yr − 1 (95th percentile).

Drivers of spatial variation in oil palm height growth across the landscape
Of the predictors included in the multiple regression model, all but TPI and aspect emerged as statistically significant in explaining variation in oil palm growth rates across the landscape (Fig. 6). Initial palm height emerged as the single strongest predictor of height growth, with smaller individuals (<4 m in height) growing an average of 1.7 m yr − 1 , while mature ones (>12 m in height) grew 1.2 m yr − 1 (Fig. 5). Canopy cover was the second strongest predictor of growth, with palms growing in denser stands exhibiting faster height growth rates. Distance from forest edges and distance to rivers were the next strongest predictors of height growth but they had opposite effects. Oil palms generally grew faster nearer forest edges and somewhat more slowly nearer watercourses. Finally, oil palms generally grew slightly faster when positioned on steeper terrain within the landscape, and at lower elevations. However, it is necessary to consider that despite being statistically significant, the effect size of these predictor variables was generally low and together they only explained around 6% of the variation in oil palm height growth rates across the SAFE landscape.

Oil palm growth during El Niño events
We found that oil palms were able to continue growing over the twoyear period between 2014 and 2016 despite experiencing drought conditions linked to El Niño. Our results are consistent with recent work showing that net primary productivity of oil palm plantations in Indonesia actually increased during the initial phase of the 2015-16 El Niño despite decreases in soil moisture availability, due to lower cloud cover and greater incoming solar radiation (Stiegler et al., 2019). It was only later, when local fires ignited during the drought period caused haze to increase, that CO 2 uptake by oil palms declined. By contrast, adjacent forests within the same landscape grew much more slowly during the El Niño period than before or after it, and even lost canopy height near edges due to tree mortality and leaf shedding (Nunes et al., 2021). Not having access to data on oil palm growth rates before and after the 2015-16 El Niño, we could not directly quantify the degree to which drought conditions impacted oil palm growth rates. However, when our results are compared to the those of Nunes et al. (2021), which quantified the response of native vegetation within the same landscape to the 2015-16 El Niño, we found that oil palm growth rates were not only much faster but also much more consistent across the landscape. This suggests that oil palm growth is limited more by light than it is water availabilitymore so than in the case of native vegetation where high vegetation density cause intense competition for water during periods of drought.  That being said, our results are at odds with the conventional wisdom that palm oil production would tend to decline during periods of drought, such as those associated with El Niño (Hsiao, 2003;Oettli et al., 2018). For example, previous work by (Woittiez et al., 2017) found that when there are water deficits >400 mm per year, oil palm yields were more than halved. Similarly, long-term observations from peninsular Malaysia, Sabah and Sarawak suggest that during extreme El Niño events like the one of 1997-98, oil palm fresh fruit bunch yields were >10% lower than in normal years (Oettli et al., 2018). What this suggests is that while oil palms may be able to maintain CO 2 uptake and vegetative growth during periods of low water availability, this does not necessarily then translate into fruit production. This highlights an important caveat to our analysis, as vegetative growth is not necessarily a direct proxy for fruit or oil yields, particularly over short-term periods characterized by anomalous climatic conditions.

Variation in oil palm growth rates across landscapes
While we anticipated that water availability would be one of the strongest predictors of palm growth across the landscape during El Niño events, oil palm growth rates actually increased with distance from rivers and decreased further away from forest edges. These trends could be due to a variety of reasons, including the fact that distance from rivers may not always represent the best proxy for soil moisture availability or as mentioned previouslythat growth may be faster when soils are not fully saturated (Stiegler et al., 2019). What our results do suggest is that common-place strategies aimed at minimising forested riparian buffer zones on oil palm plantations may be misguided. Many oil palm plantation managers view forest patches as potential sources of pests and disease and therefore would prefer to remove these where possible (Edwards et al., 2014;Salaheen and Biswas, 2019). Our results suggest that the oppositive may in fact be the case and that maintaining and expanding riparian zones could allow land to be used more efficiently for oil palm production while also preserving key habitat corridors for biodiversity and reducing run-off of soil and fertilizers into river systems. This would also have the benefit of creating a more structurally complex landscapes, with benefits for both biodiversity and people (Chaplin-Kramer et al., 2011;Edwards et al., 2014;Klein et al., 2003;Koh, 2008;Landis et al., 2003;Ricketts et al., 2004;Tscharntke et al., 2008;Woltz et al., 2012).
In terms of topographic effects on growth rates, we found that oil palms grew slightly faster at higher elevations within the SAFE landscape. A possible explanation for this is that during the hotter-thanaverage El Niño conditions, temperatures in the lowest parts of the landscape may have exceeded optimal conditions for growth. This is consistent with previous work exploring how local-scale microclimatic conditions vary across this study area, which showed that in lowland areas with low vegetation cover, high near-surface air temperatures cause atmospheric vapour pressure deficit (VPD) to exceed the point at which most tropical plant speciesincluding oil palmwill shut their stomata (Jucker et al., 2018c). By contrast, we found no relationship between growth rates and TPI, indicating that small-scale variation in terrain features such as ridges and depressions which affect water flow and retention in the soil were of little importance in determining oil palm growth. This again is in contrast to what has previously been observed for non-oil palm vegetation within this same landscape (Nunes et al., 2021;Röll et al., 2015), and suggests that atmospheric dryness linked to air temperature and VPD may be more important in determining oil palm growth than local-scale differences in soil moisture availability.
Lastly, we found that height growth was strongly linked to initial size and the cumulative vegetation cover of surrounding individuals. By far the strongest predictor of palm growth in our study was the initial size of palms, with smaller, younger palms growing significantly faster than those that were larger and already mature. This likely reflects an ontogenetic trend in growth, whereby most plant species will slow down their height growth as they approach maturity. However, previous work has shown that oil palms continue to grow almost linearly in height even by the time they reach 20-25 years of age (Pang Tan et al., 2014). An alternative reason why younger oil palms may have been able to sustain faster growth rates during the El Niño period is that they tend to have much lower water use rates compared to mature individuals (Röll et al., 2015), suggesting they would have required less water to sustain their growth. In addition to initial size, canopy cover was also a significant predictor of growth, with palms tending to grow fastest when in denser stands. This is consistent with our understanding of how plants allocate resources towards vertical growth when experiencing competition for light from neighbouring individuals. It also underscores the fact that under current planting densities of 140-160 oil palms per hectare, individual palms do not appear to be competing significantly with each other for water even during conditions of drought brought upon by El Niño.

Towards monitoring oil palm growth rates and yields using remote sensing
While our approach based on the use of repeat airborne LiDAR allowed us to successfully map the height growth of individual palms, when it came to explaining how these growth rates vary across landscapes we were only able to explain a small portion of the total variance. To some extent this likely reflects that oil palm plantations have been designed specifically to minimise this variation and standardize yields. However, it may also indicate that other underlying drivers of growth differences between individuals were not accounted for in our analysis, such as variation in fertilizer use, soil type, pest prevalence or oil palm cultivar (Manjit et al., 2014;Ruslan et al., 2019;Sundararaju and Ratnakaran, 2002;Woittiez et al., 2019).
To our knowledge our study is one of the very first to use repeat LiDAR to track the growth of individual trees at landscape scales. In terms of developing this approach further, a next step would be to test whether alternative crown delineation routines could further improve the accuracy of individual oil palm detection. A number of studies have been conducted to benchmark the accuracy of different individual tree crown delineation routines (Aubry-Kientz et al., 2019;Eysn et al., 2015), including the dalponte2016 algorithm used here (Dalponte and Coomes, 2016), AMS3D (Ferraz et al., 2016), Graph-Cut , Profiler (Hamraz et al., 2017(Hamraz et al., , 2016 and SEGMA (St-Onge et al., 2015). These comparisons have generally found dalponte2016 to be among the best performing algorithms, particularly among those which can be applied to rasterized CHM data and are therefore computationally scalable to large areas. However, beyond the comparisons we carried out as part of this study (Table S1), as far as we are aware a similar benchmarking exercise has yet to be carried out specifically for oil palms, meaning that there may be room to further improve the accuracy of the individual tree detection routine presented here. Nevertheless, any attempt to further improve the detection of individual oil palms will need to carefully weigh these improvements against the computational scalability of the approach. In particular, while point cloud based crown delineation routines like AMS3D may offer some marginal increases in the accuracy of single tree detection, they are currently computationally much too slow to run on areas spanning hundreds of square kilometres such as the SAFE project landscape considered in this analysis.
Another key area of future work lies around establishing whether there is a robust link between oil palm vegetative growth, which can be measured via remote sensing, and yields, which instead may not be possible to estimate directly from airborne imagery. This also raises the question of whether airborne LiDAR is still the best tool for this task. RGB imagery acquired using low-cost UAVs is proving increasingly promising as a tool to both delineating individual plants using deep learning (Weinstein et al., 2020) and then measuring their height using structure from motion (albeit with its limitations linked to poor terrain detection in closed canopy systems). As the resolution of these data are far higher than that typically achieved by LiDAR-based CHMs (centimetres as opposed to meters), it may even be possible to use these data to directly map oil palm fruit bunches, allowing a way to explicitly quantify both vegetative growth and yields at the same time. As our results demonstrate, there may well be opportunities for better management of oil palm plantations to achieve win-win scenarios for both agricultural production and biodiversity, particularly around the retention and expansion of riparian buffer zones. Novel remote sensing approaches allowing both growth and productivity of individual oil palms to be tracked over time can help pinpoint where these low hanging fruit lie.

Author statement
TJ designed the study. GPA and DAC provided the data. LB conducted the analysis with the assistance of TJ and MD. LB wrote the first draft of the manuscript with TJ. All authors contributed substantially to revisions.

Data availability statement
The 2014 DEM, CHM and point density data are archived online and freely available at: 10.5281/zenodo.4020697. All data needed to replicate the results of this study will be archived on Zenodo following the review of this article.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
All data used in this study are publicly available and a corresponding DOI has been added to the manuscript.