Terrestrial Laser Scanning to Detect Liana Impact on Forest Structure

Tropical forests are currently experiencing large-scale structural changes, including an increase in liana abundance and biomass. Higher liana abundance results in reduced tree growth and increased tree mortality, possibly playing an important role in the global carbon cycle. Despite the large amount of data currently available on lianas, there are not many quantitative studies on the influence of lianas on the vertical structure of the forest. We study the potential of terrestrial laser scanning (TLS) in detecting and quantifying changes in forest structure after liana cutting using a small scale removal experiment in two plots (removal plot and non-manipulated control plot) in a secondary forest in Panama. We assess the structural changes by comparing the vertical plant profiles and Canopy Height Models (CHMs) between pre-cut and post-cut scans in the removal plot. We show that TLS is able to detect the local structural changes in all the vertical strata of the plot caused by liana removal. Our study demonstrates the reproducibility of the TLS derived metrics for the same location confirming the applicability of TLS for continuous monitoring of liana removal plots to study the long-term impacts of lianas on forest structure. We therefore recommend to use TLS when implementing new large scale liana removal experiments, as the impact of lianas on forest structure will determine the aboveground competition for light between trees and lianas, which has important implications for the global carbon cycle.


Introduction
Aboveground forest structure is an important factor influencing biodiversity, net primary productivity and the carbon cycle of tropical forests.Tropical forests are undergoing large-scale structural changes due to anthropogenic disturbances such as increased atmospheric CO 2 , logging, hunting, and conversion of forested areas into agricultural lands [1][2][3].One such structural change in tropical forests is the increase in liana abundance and biomass in the Neotropics.The mechanisms that explain liana proliferation include increased atmospheric CO 2 , evaporative demands, forest fragmentation and forest disturbance [4].Lianas are woody climbing plants that use trees and other plants as structural support for ascending to the canopy.They allocate more resources to canopy development, stem and root elongation than to their structure, resulting in a high leaf to stem ratio compared to trees [5].As a result, lianas can significantly attenuate light in the forest and can contribute up to 40% of the canopy leaf cover [6].
Due to their erratic growth forms and high abudance, lianas play an important role in tropical forest ecology and management.For example, lianas with their high canopy to stem ratio can physically link trees together in the canopy, which serves as a great canopy access system for arboreal animals [7,8].Studies have shown that lianas tend to colonize tree fall gaps leading to liana-dominated patches in the forest [9].Schnitzer and Carson [10] showed that gaps maintain better liana diversity and species richness compared to shade-tolerant tree and overall plant diversity.Both of the above-mentioned factors play a vital role in logged forest.Logging a tree with inter-crown liana connections could damage the neighboring trees in the connection and thus result in large canopy gaps [11].In addition, following logging, lianas in the crown of the logged trees can rapidly regenerate and colonize the gap and the surrounding trees [12,13].
Pre-felling liana cutting has proven to be a successful management strategy to reduce the impact of liana-induced damage in a logged forest and to promote tree regeneration in the gaps [14][15][16][17].Apart from the forest management strategy, liana-cutting experiments are conducted to study their long-term impacts on tree-growth and carbon dynamics of the forest [15,18].Long-term impact of liana cutting is often measured by field observations like tree diameter increment, which is then linked to aboveground biomass using allometric equations [19].Van Der Heijden et al. [19] showed that, after three years of liana removal in large plots in Panama, the net above-ground carbon uptake of the trees was reduced by 76% in the control plots compared to the removal plots, mainly due to the liana-induced mortality of trees.Nevertheless, the changes to the vertical and horizontal structure of the forest, which is directly linked to the aboveground competition between lianas and trees, caused by liana removal is seldom studied.Pérez-Salicrup [15] reported an increase in canopy openness of up to 4%, measured by hemispherical photographs, 26 months after liana cutting.
Recent advances in remote sensing technologies have enabled us to view the forest structural complexity in new and unprecedented ways.Terrestrial laser scanning (TLS) is an active remote sensing technique and can measure various forest structural parameters with high spatial accuracy.Vertical distribution of Plant Area Index (PAI) and Plant Area Volume Density (PAVD) are key forest structural metrics directly related to the light interception, growth and primary productivity of the forest.PAI estimates derived from TLS data were found to be more robust than the PAI estimated from hemispherical photography [20,21].In addition, various studies have successfully derived vertical profiles of gap fraction, PAI and PAVD for forest ecosystems from TLS data [22][23][24].Various other metrics such as the height of maximum vegetation density [24], Canopy Height Models (CHM) [25], stem density maps [26] have also been derived from TLS data with high accuracy.
However, only a handful of studies have used TLS data for characterizing the changes in the forest structure.Calders et al. [27] used multi-temporal TLS measurements based on 48 measurement days from four sampling locations to monitor vegetation phenology by characterizing the shift in the vertical distribution of PAI.Srinivasan et al. [28] have investigated the use of multi-temporal TLS data to measure the growth of the trees by estimating the change in tree-level aboveground biomass.Liang et al. [29] detected changes to the stems over time caused by natural or other forces using bi-temporal TLS data.Olivier et al. [30] developed a method for quantifying the canopy changes, mainly the tree crown responses to gap formation by collecting TLS data before and two years after gap formation.Despite the recent technological advances and different studies highlighting the impact of lianas on the forest structure, very few studies exist that have tried to analyze the forest structural change in the 3D canopy structure caused by liana removal.A first study used a Portable Canopy LiDAR (PCL) [31] and Li-COR LAI-2000 plant canopy analyzer (LI-COR, Lincoln, NE, USA) to characterize the change in forest canopy before and each of the two years after liana cutting [6].PCL is an upward looking infrared pulsed-laser system that has to be moved through the forest in vertical transects at a constant pace and spacing.The change in the plant surface density across different vertical strata between the removal and control plots indicated the contribution of lianas to upperand mid-canopy of the forest.However, the contribution of lianas to the upper canopy of the forest was not quantified owing to the technical limitation of the PCL.Other studies have used VEGNET, a phase-based laser rangefinder that reflects the pulses at a fixed angle of 57.5 • zenith, to assess whether liana presence influences the forest successional trajectories by assessing the vertical signature of the forests of different stand age with and without lianas [32,33].
Since lianas are disturbance-adapted plants, liana abundance is likely to increase with increased forest disturbance, thereby increasing tree mortality and decreasing tree growth.Our study assesses the potential of TLS observations of liana removal experiments to better understand the role of lianas in forest functioning.Forest structural changes owing to liana proliferation could explain aboveground competition between lianas and trees, which has important implications for the forest carbon cycle [34].In our study, we use for the first time high spatial resolution 3D measurements from a multiple-return TLS to study lianas in the tropical forest.We collected bi-temporal TLS data from a removal plot before and six weeks after liana cutting and a control plot of the same size.Based on this test dataset, we evaluated the potential and limitations of high resolution TLS to observe the vertical and horizontal structure of lianas in the forest.
The specific objectives of this paper are: 1. to demonstrate the potential of TLS to detect changes in the vertical structure of the forest before and after liana removal 2. to study the reproducibility of the TLS-derived metrics by comparing the structural metrics from two time steps of the control plot.

Study Site
We conducted a liana removal experiment at Gigante Peninsula in Panama (9 • 9 N, 79 • 51 W), which is adjacent to Barro Colorado Island (BCI) and is a protected forest part of the Barro Colorado Natural Monument (BCNM).This forest is a 60-year old secondary seasonal tropical forest with recorded average annual rainfall of 2600 mm per year.The site has a distinct dry period for four months from December to April during which the rainfall seldom exceeds 100 mm/month [10,19].Gigante Peninsula has predominantly clay-rich, highly weathered Oxisols.Tree density and liana density of the site were found to be 3600 stems and 2000 stems ≥1 cm per ha in 2008 [18].The map of the study site with the control and removal plots marked as red dots is shown in Figure 1.

Experimental Set-Up
We set up two 10 m × 10 m plots, one removal plot and one control plot in February 2016.We measured all trees with diameter ≥10 cm at breast height (1.3 m from ground) and all lianas with diameter ≥1 cm following the liana census protocol in [35] before liana removal.The removal and control plots were approximately 1.5 km far from each other.Since cutting down lianas could release the neighboring trees from below-ground competition, we chose the control plot far from the removal plot [34].The removal plot had a total of five trees and 66 lianas that were rooted in the plot.Eighteen lianas that were rooted in the plot ascended into the canopy using host trees outside of the 10 m × 10 m plot.We recorded the liana load on trees by counting the number of liana stems that rooted within the plot and attached themselves to the tree trunk or branches in the removal plot.It is possible that some lianas root outside of the plot and go into the canopy of the trees inside the plot.However, it is practically very difficult to follow these lianas as they can root far from the host tree.Hence, only the lianas that root inside the plot are included in the inventory as mentioned in [35].The control plot had a total of four trees and 45 lianas rooted in the plot.In Table 1, we provide the list of tree diameter, height and liana load of each tree in the removal plot and the tree diameter and height from the control plot.We derived the height of the trees using the region-specific diameter at breast height (DBH) to height allometric equation listed in the study of Van Der Heijden et al. [19].The height derived from allometric equations do not represent the true height of the trees but rather give an estimate of the tree height.Trees 2 and 3 were actually the same individual that branched below 1 m from the ground.As a result, we counted the lianas load on one tree also for the other, as it was difficult to distinguish the canopy of these two trees from the ground.Lianas were cut near the forest floor using a machete and most of the liana stems in the understorey were moved without damaging the tree crowns on 23 February 2016 as it can be seen in Figure 2. Lianas in the area of 2 m around the plot were also cut to avoid the edge effects.

LiDAR Data
We collected TLS data from the two 10 m × 10 m plots, using four scan locations at the corners of the plot.TLS data was collected twice in the removal plot, once on 22 February 2016 one day before liana removal and again on 3 April 2016, six weeks after liana removal, assuming that liana leaves would have mostly fallen off the canopy in the six week time period [6].We scanned the control plot of size 10 m × 10 m twice with the same time gap of six weeks (26 February 2016-4 April 2016) to account for natural leaf abscission.
We used a RIEGL VZ400 terrestrial laser scanner (RIEGL Laser Measurement Systems GmbH, Horn, Austria), which is a multiple return time-of-flight based scanner using a narrow infrared laser beam of wavelength 1550 nm and a beam divergence of 0.35 mrad.We mounted the scanner on a tripod at approximately 1.3 m from the ground and used an angular sampling resolution of 0.02 • .The acquisition time for one full scan with an angular resolution of 0.02 • is 12 min 19.6 s. Figure 3 illustrates the view from one scanning position from the removal plot before and after liana removal in a cylindrical projection.Figure 3b   As the scanner has a zenith range of only 100 • from 30 • to 130 • , to get a full hemispherical view of the canopy, the scanner is tilted to 90 • from the vertical axis for an additional tilted scan.We distributed reflective targets in the field to register the upright scan with the corresponding tilt scan and also to co-register all four scan locations.The registration of all four positions for the removal and control plots was done using the RISCAN Pro software (version 2.5.3,RIEGL Laser Measurement Systems GmbH, Horn, Austria) provided by Riegl.After registration, the raw TLS data were exported as ASCII files for deriving CHMs and stem maps.In addition, the raw TLS data were converted into an open source Sorted Pulse Data (SPD) format and analysed using the open source LiDAR data processing library PyLidar for deriving the vertical profiles of Pgap, PAI and PAVD [36,37].

Co-Registration of Bi-Temporal Data
We co-registered the TLS data before and after liana removal in the removal plot and two scans from the control plot to facilitate the comparison of the all TLS-derived metrics from two different time steps.The point cloud of the region of interest from the two scans of removal and control plot were exported as an ASCII file from RiSCAN pro.The co-registration was done on these ASCII files using an Iterative Closest Point (ICP) algorithm implemented in CloudCompare (version 2.8.1, CloudCompare, GPL software) [38].The detailed steps of co-registration of the point clouds are described in Appendix A.2.We refer to the normalized z-value as height in this article.Since the topography of the plots studied are relatively flat, we went for a simple plane-fitting normalization as explained in [24].

Vertical Plant Profiles
The frequency distribution of all the returns with respect to height for the bi-temporal data of the removal and the control plot were estimated using a Gaussian kernel density estimator with a bandwidth of 0.5.We also estimated the frequency distribution of the different returns (first, second, third and fourth) with respect to height for the pre-and post-cut data of the removal plot using the kernel density estimator with the same bandwidth of 0.5.
We performed a two-sided Kolmogorov-Smirnov test (KS test) to test the equality of the relative frequency distribution of height of the two time steps for the removal as well as control plot [39].We chose this test as it does not assume that the data are sampled from a defined distribution like Gaussian distribution.
We described the vertical plant profiles using the vertically resolved gap probability, PAI and PAVD derived from the TLS data.As mentioned in the Section 2.3, we used the open source python library PyLidar for deriving the vertical profiles.We followed the method described in Calders et al. [24] to estimate the Pgap, PAI and PAVD profiles.PAI(z), PAI as a function of height z, is calculated through the estimates of Pgap at the corresponding height.PAVD(z) is calculated based on vertically averaged PAI at 57.5 • zenith and Pgap(z).The basic steps and formulas are described in Appendix A.1.For the removal and control plot, we derived the profiles for the specific azimuth region of one scan location shown in Figures 3 and 4, respectively, that fell within the 10 m × 10 m plot.We did not derive the profiles for the data from other scan locations as the scans had many objects too close to the scanner (<1.5 m from the scanner), especially in pre-cut scan, leading to artificial gaps resulting in underestimated PAI.

Nearest Neighbor Distance
The angular resolution (0.02 • in this case) and the range, the distance between the target and the instrument, determine the vertically resolved Average Nearest Neighbor (ANN) distance for the point cloud.The minimum of the average nearest neighbor distance at a particular range is the smallest distance resolvable by the TLS instrument (d min ), which depends on instrument-specific beam divergence [40].Vertically resolved ANN is calculated as in Equation (1): NN is the nearest neighbor distance for a point in the point cloud.The distance between two points in the point cloud for a plot is defined by the Euclidean distance between the two points in 3D.N(z) is the total number of points in the height interval between z and z + dz.For every 1 m z-bin, we calculate the average distance to four nearest neighbors for all the points falling in the bin.
The ANN(z) values indicate the average distance between the topologically connected points from forest floor to canopy.In forests with dense vegetation, ANN(z) is not only influenced by d min , but also by the spatial distribution of the forest elements occluding the other elements with increasing height [40].In our study, we compare the ANN(z) between the pre-and post-liana cut TLS data from the removal plot to assess the complex spatial distribution of liana stems and leaves in the plot.We also compare the ANN(z) between the multiple return and the first return TLS data to assess the advantage of using multiple-return scanners over the first-return scanners for studying lianas in the tropics.

Canopy Height Models
Canopy Height Models describe the horizontal distribution of the height of a forest canopy as a three-dimensional surface [41].In this study, we generated CHMs to compare the change in height in the observed canopy of the co-registered point clouds in the removal plot before and after liana removal.We also generated CHMs for the two point clouds from the control plot and compare it to that of the removal plot.Bi-temporal TLS data from the removal and control plots were co-registered as described in Appendix A.2.We derived CHMs by selecting the points of highest z-value within each 50 cm × 50 cm x,y grid.After generating CHMs for the two point clouds, we calculated the absolute difference between the corresponding grids of the bi-temporal CHMs to compare the change in the spatial distribution of the observed canopy height for both the removal plot and control plot.

Vertical Plant Profiles
The frequency distribution of returns with respect to height shows major differences in the number of returns from understorey and the upper canopy of the bi-temporal scans in the removal plot (Figure 5a).After cutting lianas, we found a 10% decrease in the absolute number of returns from the first five meters and a 56% increase in the absolute number of returns above 15 m post liana cutting.The mid-canopy shows very little change before and after liana cutting in terms of the vertical distribution of returns.This is likely because the liana stems in the mid-canopy remained, whereas the liana stems in the understorey were removed while cutting and the liana leaves in the upper canopy had mostly fallen six weeks after liana cutting.In the pre-cut data, there were fewer returns from the upper canopy with lower canopy height compared to the post-cut data because of the occlusion of liana leaves obstructing the pulses emitted.
The pre-and post-cut frequency distributions of height were significantly different in the removal plot (KS statistic = 0.16, P < 0.001), but not in the control plot as it can be seen from Figure 5b (KS statistic = 0.05, P < 0.144).
We found differences in the frequency distribution of height of the different returns in the understorey and the canopy between the pre-and post-cut TLS data (Figure 6).For instance, after liana cutting, the absolute number of returns decreased by 60% in the understorey up to five meters and increased by 70% above 15 m in the post-cut TLS data when only the last returns (fourth return) were considered.Results of the frequency distribution of heights of different returns for the bi-temporal data of the control plot show similar results as in Figure 5b and are therefore not discussed further.However, the results are included in Appendix A.3 for reference.We derived the vertical profiles of PAVD with both first-return and multiple-return TLS data for the pre-and post-cut scans of the removal plot (Figure 7).The profiles were generated for the azimuth region shown in Figure 3.The total PAI decreased from 6.5 to 6.3 after liana cutting.We found 48% decrease in the first five meters, 18% decrease at 15 m and 90% increase above 16 m in PAVD post-liana cutting when using multiple returns.Both the first return and multiple return PAVD profiles show very little change in the mid-canopy before and after liana cutting similar to the vertical distribution of returns shown in Figure 5a.Both the profiles show increased PAVD beyond 16 m.This increase in PAVD beyond 16 m mainly corresponds to the occlusion of liana leaves in the canopy at about 15 m.Our results demonstrate that liana significantly increases occlusion in the upper canopy owing to their high leaf to stem ratio [5].
Similar to Calders et al. [24], we observed a larger number of laser pulse returns at greater canopy height when using multiple returns as shown in Figures 6 and 7.This highlights the advantage of using a multiple return instrument in forest ecosystems with dense understory and high liana abundance.We recorded 270% more returns from the canopy (≥15 m) when using multiple returns compared to first returns in the pre-cut and post-cut TLS data.In addition, the highest return of pulses for the pre-cut data went down to 18 m from 20 m when only first returns were considered.For the post-cut data, the highest return went down to 19.5 from 22 m.A two meter increase in TLS estimated height post liana cutting clearly indicates that the complex, spatial distribution of different forest elements results in increased occlusion effects with increased height.The height of the largest tree according to region-specific diameter to height allometric equation was 25.19 m and 26.11 m in the removal and control plot, respectively.The maximum height reached by the multiple return TLS used in the study was 22 m (after liana cutting) and 24 m in the removal and control plot, respectively.There are two plausible explanations for this height difference of 3 m in the removal plot and 2 m in the control plot.One possible explanation is that the TLS is underestimating the maximum canopy height owing to tree top occlusions.Tree height measurements from vertex rangefinders have proved to be more accurate in some forest ecosystems compared to TLS and other traditional height measurement methods [42][43][44].However, this is not always the case as TLS has been proved to reach the top of the canopy easily in other forest ecosystems, where the TLS-derived tree height was compared to true tree height from destructive measurements [45].The other explanation is that the allometric equation used is overestimating the tree height.The PAVD profile with only first returns demonstrates more vegetation closer to the scanner in the understory than in the upper canopy compared to the multiple-return PAVD profile.The peak PAVD height before liana cutting remains unchanged at 15.5 m for both first and multiple return profiles.This demonstrates that both the first return and multiple return profiles agree on the height of higher concentration of liana leaves.However, the peak PAVD height after liana cutting went down to 13.5 m from 15 m demonstrating the downside of first return TLS instruments of having a low amount of returns from the canopy despite the mitigated occlusion effects from the loss of other canopy elements.
The bi-temporal vertical profiles of PAVD of the control plot are shown in Figure 8a. Figure 8b shows the box plot of absolute difference between the PAVD of two time steps of the control and removal plot.As it can be seen from Figure 8b, the bi-temporal PAVD profiles of the control plot did not show any major shifts.Our results demonstrate not only the reproducibility of the TLS-derived metrics to do repeated measurements to study the temporal dynamics of tropical forest structure as a function of height, but also the potential of TLS to determine the changes in the vertical profile followed by liana removal.Reproducibility of TLS-derived metrics has also been confirmed in a broadleaf deciduous forest by a spring phenology monitoring study of Calders et al. [27].
We only used one location to derive the vertical profiles of PAVD in the removal plot and control plot.This was mainly because of many objects being close to the scanner (<1.5 m from scanner) especially in the pre-cut scan indicated by the blue zones in Figure 9. Including the blue zones in the zenith band of 55-60 • for the estimation of Pgap will result in an overestimated Pgap and thus an underestimated PAI.One way to overcome this is to scan from a location with no interfering vegetation too close to the scanner or to enable the near range activation mode (which is scanner specific, 0.5 m for our instrument).The other way is to correct for these effects using voxel-based methods [46,47], which will be an area of future study.

Nearest Neighbor Distance
Comparison of the multiple-return and the first return ANN (z) for the pre-and post-cut TLS data of the removal plot and the bi-temporal TLS data of the control plot are shown in Figure 10.Lower ANN value for a particular height compared to the corresponding d min indicates theoretical oversampling and higher ANN value indicates theoretical undersampling.The line representing the laser beam diameter at different heights shows what theoretically can be resolved by the laser beam.The higher the deviation from this line, the higher is the occlusion and lower is the quality of the data.
Pre-cut TLS data clearly shows higher effects of occlusion than the post-cut TLS data in the upper canopy, especially after 15 m.These effects are intensified when using only first returns as the ANN increases by 2 cm from 1 cm and 6 mm from 4 mm between 15 m and 17 m in the pre-cut data and post-cut data, respectively.It is evident from PAVD profiles of Figure 7 that the highest concentration of liana leaves is located at about 15 m leading to an increased ANN above 15 m in the pre-cut data.This implies that, in the first return pre-cut data, the size of the smallest object that could be resolved at 17 m is at least 3 cm and there are not enough topologically connected objects beyond 17 m, whereas with multiple returns the minimum size of the object resolvable at 18 m is 2.5 cm.As this trend is only expected to increase when we move higher up in the canopy, multiple-return TLS has an advantage over first-return TLS to study impact of lianas on forest structure as liana leaves are mainly in the canopy.
Bi-temporal TLS data from the control plot shows no significant difference in the vertical ANN profile.The multiple return ANN at 21 m remains at 4.1 cm and 4.2 cm for the TLS data collected on 26 February 2016 and 4 April 2016, respectively.The same is true for the first return ANN at 19 m, which stays at 3 cm and 2.9 cm for the TLS data on 26 February 2016 and 4 April 2016, respectively.This again demonstrates the potential of TLS to produce reproducible results and to determine the changes in the vertical profile of the forest.However, this also shows the occlusion effects of dense vegetation beyond 15 m, which should be addressed in the upcoming studies focusing on lianas in the tropical biomes.

Canopy Height Models
The difference between the pre-and post-cut CHMs of the removal plot and the bi-temporal scans' CHMs of the control plot are illustrated in Figure 11.We observed an overall increase in canopy height in the CHM after liana cutting with 82% of the grids seeing an increase in height between 0 to 6 m and 18% of the grids losing vegetation between 0 and 4 m.There were some local changes in the CHM between two time steps in the control plot.The observed increase in canopy height was approximately the same as the observed decrease in the canopy height in the control plot with 54% of the grids of the TLS data from 4 April 2016 gaining height between 0 and 4 m and 46% of the grids losing height between 0 and 4 m.The average observed increase in canopy height was 1.5 m in the removal plot and 0.4 m in the control plot, whereas the average observed height loss was 0.5 m for both the removal and control plot.
It can be seen from Figure 11 that contiguous grids with observed increase in canopy height ranging between 1 and 6 m are co-located with the location of trees 4 and 5 in the removal plot.We interpreted the virtual height gain in these contiguous grids as loss of liana leaves in those grids.Clark and Clark [48] estimated the liana load on trees by two measures: the number of liana stems attached to the tree trunk and Crown Occupancy Index (COI).The COI is estimated on a scale of five based on the percentage of the tree crown occupied by liana leaves (0, 1-25%, 26-50%, 51-75%, 76-100%).They showed that the COI of trees were positively correlated with the basal area of liana load attached to their trunks.We found evidence for their hypothesis as trees 4 and 5 had higher basal area of liana load than the other trees in the plot corresponding to loss of more liana leaves in the canopy.Figure 12 shows the height at which these changes occurred pre-and post-liana cutting in the removal plot.We showed that most of the observed height gain of more than 2 m occur beyond 14 m in the pre-cut data, which is also evident from PAVD profiles in Figure 7.As tree 5 is a small tree with 15.3 cm DBH and 14.9 m height, it is reasonable to assume that most of the height gain occurred within the crown of tree 4, the largest tree of the plot.Our results are in agreement with the other studies, which showed lianas compete more with upper canopy trees for light by colonizing the tree crowns [49,50] and the study of Clark and Clark [48] that showed the liana loads to be positively correlated with tree diameter.
Difference between CHMs of the control plot showed local changes in the canopy.These changes might be either due to the displacement of vegetation or due to loss of some leaves or branches within the six-week period as indicated by Figure 4.However, as Figure 11b,d indicate, the scatterplots confirm the height agreement of bi-temporal data of the control plot as opposed to the bi-temporal data of the removal plot.The TLS instrument used in the study shows the potential for capturing the spatial distribution of changes in the canopy in high detail, which has important implications for the long-term monitoring of liana removal experiments with both removal and control plots [6].

Limitations of the Study
As explained in the previous sections, occlusion is a challenging problem when deploying TLS in complex environments like tropical forests.The observed increase in canopy height after liana cutting in Figure 11a does not correspond to the real increase in canopy height but rather indicates the occlusion effects of liana leaves in the pre-liana cut CHM.We collected TLS data only from the four corners of the plot.One of the ways to overcome the occlusion is to have a better sampling strategy by increasing the number of scanning locations.Increasing the number of scans around the plot enables us to cover the plot from multiple viewing angles.In addition, scanning the plot from longer distances (approximately 5-10 m from the plot edge) could result in getting better views of the canopy [51].Our study was able to detect the changes in the vertical structure of the forest, but a better sampling strategy could result in data with better quality enabling the quantification of the changes with higher accuracy.In addition, future development of TLS instruments with smaller beam divergence could help in resolving smaller objects at greater heights in the canopy, thereby mitigating the occlusion effects.
Another important limitation of the study is the plot size and sample size.We used a small sample size because our study is a first exploratory effort to test the potential of TLS to study the effect of lianas on forest structure.We tried to balance the number of lianas that we sacrificed for this study by cutting enough lianas to rigorously see the effect of liana cutting on forest structure, but not so many that it would have a huge impact on natural forest regeneration and functioning.Though our study cannot conclude about the ecological impact of lianas on forest structure, the changes observed in the PAVD profiles, ANN(z) profiles and the CHMs of the removal plot derived from the TLS data have important implications for the future studies of liana impact on the forest structure.Collecting high quality TLS data for large areas in dense forest ecosystems can be time-consuming.For instance, scanning a 1 ha plot with scan positions located every 10 m and few scan positions located outside of the plot boundaries in a dense tropical forest could take up to five days including the time taken for preparing the plot (distributing the control reflective targets for registration).However, such high detailed data can yield answers to some of the unsolved ecological questions about liana influence on forest structure.

Conclusions
As lianas strongly compete with trees for light, studying the changes in the vertical forest structure with and without lianas over the long term could explain how lianas impact the forest structure.We demonstrate the potential of TLS to study the changes in all the vertical strata of the forest caused by liana removal.Our results demonstrate the advantages of having multiple returns over first returns to quantify changes in the upper canopy of the forest, especially in environments where the understorey is dense or lianas are abundant.Our results also highlight the reproducibility of TLS derived metrics to do repeated long-term measurements of the liana removal and control plots.Though collecting TLS data from large areas in complex environments can be time-consuming, lessons learned from the current small scale study will enable us to collect high quality data in future.We highly recommend to include TLS observations in future large scale liana removal experiments in the long term to quantitatively study the impact of lianas on the forest structure.
4. We then applied the transformation matrix from the manual registration to the whole point cloud to be aligned.5. We used the ICP algorithm [53] implemented in CloudCompare for fine registration after the first coarse manual registration.We selected all the points from the ground up to 4 m for fine registration.6.We applied the transformation matrix that resulted from the ICP fine registration to the whole point cloud to be aligned.
Co-registration of the bi-temporal TLS data was achieved with an RMSE of 0.08 m in the removal plot and 0.06 m in the control plot, respectively.Figure A2 shows the result of co-registration using a tree trunk from the removal plot.

Appendix A.3. Additional Results
Result of the frequency distribution of height of different returns for the bi-temporal data of the control plot, which shows similar results as in Figure 5b, is shown in Figure A3.

Figure 1 .
Figure 1.Map of the study site with the location of the removal and control plots marked with the red dot.The study site is located in Gigante Peninsula in Panama (9 • 9 N, 79 • 51 W) adjacent to Barro Colorado Island (BCI).

Figure 2 .
Figure 2. Comparison of the removal plot before (left) and after (right) liana removal.

Figure 3 .
Figure 3.View of terrestrial LiDAR data from one scan position of the removal plot with azimuth ranging from 120 • to 240 • and zenith ranging from 30 • to 130 • .The data is colored according to the range of the laser pulse.

Figure 4 .
Figure 4. View of terrestrial LiDAR data from one scan position of the control plot with azimuth ranging from 0 • to 100 • and zenith ranging from 30 • to 130 • .The data is colored according to the range of the laser pulse.

Figure 5 .
Figure 5.Comparison of frequency distribution with Gaussian kernel density estimation of all vertical returns combined for the bi-temporal TLS data from the removal (a) and control plot (b).

Figure 6 .
Figure 6.Comparison of the frequency distribution with Gaussian kernel density estimation of different vertical returns with height between the pre-cut (a) and post-cut (b) data of the removal plot.

Figure 7 .
Figure 7.Comparison of the first-return and multiple-return Plant Area Volume Density (PAVD) profiles of the removal plot.

Figure 8 .
Figure 8.(a) shows the bi-temporal Plant Area Volume Density (PAVD) profiles of the control plot and (b) shows the absolute PAVD difference box plot of the control and removal plot, with t 2 being the data from second time step and t 1 , the first time step.

Figure 9 .
Figure 9. TLS data before liana removal with azimuth ranging from 0 • to 360 • and zenith ranging from 30 • to 130 • .The data is colored according to range of the laser pulse with black being low and white being high.The areas with blue indicates regions of no returns either because the object was too close to the scanner (<1.5 m) or because of a gap.

Figure 10 .
Figure 10.Comparison of multiple return vs. first return vertically resolved average nearest neighbor distance (ANN(z)) for pre-and post-cut TLS data of the removal plot (a) and bi-temporal data of the control plot (b).(• Laser beam diameter, First return and + Multiple return).

Figure 11 .
Figure 11.Difference map of the Canopy Height Models (CHM) with circles representing the location of trees ≥10 cm in DBH.(a) removal plot and (c) control plot.The size of the circles correspond to the size of the trees.The color scale represents the change in height in meters between the two time steps.Scatter plot showing the CHM difference between bi-temporal TLS data for removal plot (b) and control plot (d).

Figure 12 .
Figure 12.Differences in Canopy Height Models (CHM) for CHM post-cut − CHM pre-cut ≥ 2 m in the removal plot.

Figure A2 .
Figure A2.Illustration of co-registration using a tree trunk from the removal plot.Black: Pre-cut data and Red: Post-cut data.

Figure A3 .
Figure A3.Comparison of the frequency distribution with Gaussian kernel density estimation of different vertical returns with height between 26 February 2016 (a) and 4 April 2016 (b) data of the control plot.

Table 1 .
Overview of trees, diameter and number of liana stems on each trees in the removal and control plot.Trees 2 and 3 are the same individual that branched below 1 m from the ground.