Urban Forest Growth and Gap Dynamics Detected by Yearly Repeated Airborne Light Detection and Ranging (LiDAR): A Case Study of Cheonan, South Korea

: Understanding forest dynamics is important for assessing the health of urban forests, which experience various disturbances, both natural (e.g., treefall events) and artiﬁcial (e.g., making space for agricultural ﬁelds). Therefore, quantifying three-dimensional changes in canopies is a helpful way to manage and understand urban forests better. Multitemporal airborne light detection and ranging (LiDAR) datasets enable us to quantify the vertical and lateral growth of trees across a landscape scale. The goal of this study is to assess the annual changes in the 3-D structures of canopies and forest gaps in an urban forest using annual airborne LiDAR datasets for 2012–2015. The canopies were classiﬁed as high canopies and low canopies by a 5 m height threshold. Then, we generated pixel- and plot-level canopy height models and conducted change detection annually. The vertical growth rates and leaf area index showed consistent values year by year in both canopies, while the spatial distributions of the canopy and leaf area proﬁle (e.g., leaf area density) showed inconsistent changes each year in both canopies. In total, high canopies expanded their foliage from 12 m height, while forest gap edge canopies (including low canopies) expanded their canopies from 5 m height. Annual change detection with LiDAR datasets might inform about both steady growth rates and di ﬀ erent characteristics in the changes of vertical canopy structures for both high and low canopies in urban forests.


Introduction
Quantifying forest gaps is essential for monitoring the stability of the forest structure because these disturbances can change the light environment and drive forest dynamics [1][2][3][4]. Physical disturbances in forest canopies usually result from deforestation (e.g., making space for agricultural fields) and tree-fall events [1]. Forest gaps could be characterized by their size distribution. The frequency of forest gap areas in natural forests follows a power law distribution. Usually, small openings (formed by wind damage or tree mortality) in natural forests dominate areas in size distribution and show larger scaling exponents, while frequent large openings show smaller scaling exponents. Therefore, by comparing the scaling exponents with other forests, it is possible to compare the degrees of disturbances [2,3].
Urban forests are usually fragile to disturbances because the dense human population leads to severe conditions that can damage the forest directly (e.g., land-cover changes from forested to Thus, we estimated three-dimensional canopy changes in an urban forest using four years of annual LiDAR datasets. To assess how much urban forest structures have changed, we estimated: 1) the distribution of the growth area and damaged areas (e.g., vertical and lateral growth as well as damaged areas), 2) changes in the vertical leaf area density profile, and 3) the dynamics of opening and closing forest canopies during 2012-2015. Accordingly, this study addressed two questions: 1) What are the differences between canopy structure changes derived from annual change detections and three-year interval change detection? 2) What are the characteristics of structural changes by different canopy classes (e.g., high canopies and lower canopies) in urban forests?

Study Site
The study was conducted at Mt. Bongseo (36.82°N, 127.12°E, altitude: 158 m, area: 124 ha), a remnant patch after urban development located in the center of metropolitan Cheonan, Chungcheongnam-Do, Republic of Korea. The study site is covered with a mixed forest dominated by Quercus acutissima (percent area: 44.8%), Pinus rigida (29.7%), and Larix leptolepis (7.4%). The average age class (as designated by the Korea Forest Service, 2015, http://www.forest.go.kr/) of the forest is 3.8 (i.e., 30-40 years old; the third class). Both managed and unmanaged forest areas can be found within the study site. The forest is easily accessed by citizens, and local amenities such as health facilities are present within the forest. Moreover, some of the land in the site is used as a cemetery and some have agricultural activity, which is private land (Figure 1).

Field Survey
To classify open canopy types (e.g., natural treefall or thinning, private graveyards, facilities, and agricultural activities), we conducted field trips on 7 October 2016, 26 August 2017, and 26

Field Survey
To classify open canopy types (e.g., natural treefall or thinning, private graveyards, facilities, and agricultural activities), we conducted field trips on 7 October 2016, 26 August 2017, and 26 September 2017. Moreover, to determine how much canopy should be classified as open canopy, we conducted simultaneous localizing and mapping (SLAM) mobile LiDAR (Stencil, Kaarta, PA, USA) scanning under the tree-fall area ( Figure 2) at the site on 26 March 2018, when the trees were defoliated. We generated 0.25 × 0.25 m 2 digital surface models using the maximum height (HMAX) of all mobile LiDAR datasets with CloudCompare software (https://www.danielgm.net/cc/) and generated canopy height models (CHMs) by subtracting the digital terrain models from the digital surface models. Finally, we calculated the gap size in ArcGIS ver. 10.2.2 (ESRI, Redlands, CA, USA). The average area of the 19 open canopies in the nine sampled plots was calculated as 12.48 ± 9.76 m 2 with consideration of the foliated season and previous research by Yamamoto [43] and Runkle [44], who did not consider very small openings (<5 m 2 ) as gaps. For the convenience of calculation, we classified gaps as locations where the open area was 10 m 2 or more ( Figure 3).
September 2017. Moreover, to determine how much canopy should be classified as open canopy, we conducted simultaneous localizing and mapping (SLAM) mobile LiDAR (Stencil, Kaarta, PA, USA) scanning under the tree-fall area ( Figure 2) at the site on 26 March 2018, when the trees were defoliated. We generated 0.25 × 0.25 m 2 digital surface models using the maximum height (HMAX) of all mobile LiDAR datasets with CloudCompare software (https://www.danielgm.net/cc/) and generated canopy height models (CHMs) by subtracting the digital terrain models from the digital surface models. Finally, we calculated the gap size in ArcGIS ver. 10.2.2 (ESRI, Redlands, CA, USA). The average area of the 19 open canopies in the nine sampled plots was calculated as 12.48 ± 9.76 m 2 with consideration of the foliated season and previous research by Yamamoto [43] and Runkle [44], who did not consider very small openings (<5 m 2 ) as gaps. For the convenience of calculation, we classified gaps as locations where the open area was 10 m 2 or more ( Figure 3).

Airborne Light Detection and Ranging (LiDAR) Dataset Acquisition and Registration
The airborne LiDAR datasets were acquired on 9 October 2012 (day of the year (DOY): 283), 16 May 2013 (DOY: 136), 16 June 2014 (DOY: 167), and 30 October 2015 (DOY: 303) from IGI LiteMapper 6800 (Samah Aerial Survey), which were all during the foliated season. LiDAR data were acquired in different months. However, Song and Ryu [42] reported relatively stable seasonal leaf area index (LAI) values in a forest (in the same climate zone as our study site) in South Korea. Our study site experienced the maximum number of leaves that turned red in the first and last study year on 28 October 2012 and 26 October 2015, respectively (Korea Meteorological Administration, www.kma.go.kr). Therefore, we assumed that our data were acquired in the full-leaved seasons.
The study site was surveyed following eight trajectories ( Figure 1) by overlapping 50% of the trajectory edges at an altitude of 1000 m. The beam divergence was 0.3 mrad, and field of view was 60°. All annual point densities were higher than 8 points/m 2 , which were sufficiently dense and similar in order to effectively compare datasets. The airborne LiDAR datasets were preprocessed using algorithms built into the TerraScan (Terrasolid, Helsinki, Finland) software on the MicroStation (Bentley) platform. First, we filtered noise at low points and aerial points. Second, we used classification algorithms to classify the points into ground and nonground classes. We minimized measurement errors by selecting 190 planar areas over the entire year in the digital surface model generated from the point cloud datasets for 2012, 2013, 2014, and 2015 and matched the heights with the average elevation for 2012. Finally, the height was adjusted based on the ground control points (triangulation point: 36°49′02.09″, 127°07′27.31″, ellipsoidal height: 182.17 m; unified control point: 36°48′50.29″, 127°06′55.92″, ellipsoidal height: 906,045 m; National Geographic Information Institute) and points on the roof of a building detected for all years. We selected 70 new random 16 m 2 plots, including building rooftops and roads where the heights were consistently adjacent to the site, to compare the 2012-2014 airborne LiDAR datasets with the 2015 airborne LiDAR dataset. The average differences from the 2015 dataset were -0.01 (2012), -0.04 (2013), and -0.05 m (2014). Although the validation results showed that the heights in 2015 were slightly greater than those in the other years, we concluded that accuracy anomalies in the dataset within 0.1 m were acceptable.

Airborne Light Detection and Ranging (LiDAR) Dataset Acquisition and Registration
The airborne LiDAR datasets were acquired , which were all during the foliated season. LiDAR data were acquired in different months. However, Song and Ryu [42] reported relatively stable seasonal leaf area index (LAI) values in a forest (in the same climate zone as our study site) in South Korea. Our study site experienced the maximum number of leaves that turned red in the first and last study year on 28 October 2012 and 26 October 2015, respectively (Korea Meteorological Administration, www.kma.go.kr). Therefore, we assumed that our data were acquired in the full-leaved seasons.
The study site was surveyed following eight trajectories ( Figure 1) by overlapping 50% of the trajectory edges at an altitude of 1000 m. The beam divergence was 0.3 mrad, and field of view was 60 • . All annual point densities were higher than 8 points/m 2 , which were sufficiently dense and similar in order to effectively compare datasets. The airborne LiDAR datasets were preprocessed using algorithms built into the TerraScan (Terrasolid, Helsinki, Finland) software on the MicroStation (Bentley) platform. First, we filtered noise at low points and aerial points. Second, we used classification algorithms to classify the points into ground and nonground classes.

Generation of Height Models and Change Detection
Since the research site is mountainous, and sloped areas may distort the vertical and horizontal locations of the tree canopy and apex [45,46], we did not consider segmenting the tree apices. Therefore, we simply calculated the vertical differences in CHMs for 2012-2015 that were generated from the LiDAR datasets and estimated the changes in the canopy throughout the study period. Since quantifying the changes of canopy opening areas and identifying the growth direction (i.e., vertical or lateral direction) both require high-resolution CHMs (e.g., [24,27,28]), and the beam divergence of the LiDAR sensor was 0.3 mrad, which is a 0.3 m footprint at the ground level (flight altitude: 1000 m), we assumed that 0.25 × 0.25 m 2 grid resolution would be adequate. After preprocessing the LiDAR datasets into 0.25 × 0.25 m 2 (grid level survey) digital surface models and digital terrain models, CHMs were generated using the TerraScan software ( Figure 4a). The annual average canopy heights for 2012-2015 were 11.09, 11.02, 12.29, and 12.12 m, respectively. The CHMs for each year were classified into two classes: high canopy (HC), with canopy height >5 m, and low canopy (LC), with canopy height ≤5 m ( Figure 3). In addition, we included the canopy closure (CaC) subclass, which represents areas with LCs in earlier years that became HCs in later years. Based on the estimated growth rates outlined by Song et al. [24], we assumed that the average annual growth rate of mature trees could not exceed 50 cm/year. Therefore, after calculating the changes between the two years, we classified the growth regions into four subclasses: vertical growth (VGr) areas with annual mean growth rate ≤50 cm, lateral growth (LGr) areas with annual mean growth rate >50 cm; damaged areas with negative growth rates, and no-change areas with -10 to 10 cm growth rates.
Since grid size can affect the accuracy of HMAX estimates [47], and canopy metrics based on vertical profiles can be used to accurately estimate forest structure [34,48], we generated a larger size of grids for the plot level survey using hexagonal grids with an area of 16.24 m 2 and edge length of 2.5 m (approximately the mean size of a closed gap area during the study periods) to estimate the changes in vertical profiles, leaf area density (LAD), and canopy complexity.

Generation of Height Models and Change Detection
Since the research site is mountainous, and sloped areas may distort the vertical and horizontal locations of the tree canopy and apex [45,46], we did not consider segmenting the tree apices. Therefore, we simply calculated the vertical differences in CHMs for 2012-2015 that were generated from the LiDAR datasets and estimated the changes in the canopy throughout the study period. Since quantifying the changes of canopy opening areas and identifying the growth direction (i.e., vertical or lateral direction) both require high-resolution CHMs (e.g., [24,27,28]), and the beam divergence of the LiDAR sensor was 0.3 mrad, which is a 0.3 m footprint at the ground level (flight altitude: 1000 m), we assumed that 0.25 × 0.25 m 2 grid resolution would be adequate. After preprocessing the LiDAR datasets into 0.25 × 0.25 m 2 (grid level survey) digital surface models and digital terrain models, CHMs were generated using the TerraScan software ( Figure 4a). The annual average canopy heights for 2012-2015 were 11.09, 11.02, 12.29, and 12.12 m, respectively. The CHMs for each year were classified into two classes: high canopy (HC), with canopy height >5 m, and low canopy (LC), with canopy height ≤5 m ( Figure 3). In addition, we included the canopy closure (CaC) subclass, which represents areas with LCs in earlier years that became HCs in later years. Based on the estimated growth rates outlined by Song et al. [24], we assumed that the average annual growth rate of mature trees could not exceed 50 cm/year. Therefore, after calculating the changes between the two years, we classified the growth regions into four subclasses: vertical growth (VGr) areas with annual mean growth rate ≤50 cm, lateral growth (LGr) areas with annual mean growth rate >50 cm; damaged areas with negative growth rates, and no-change areas with -10 to 10 cm growth rates.
Since grid size can affect the accuracy of HMAX estimates [47], and canopy metrics based on vertical profiles can be used to accurately estimate forest structure [34,48], we generated a larger size of grids for the plot level survey using hexagonal grids with an area of 16.24 m 2 and edge length of 2.5 m (approximately the mean size of a closed gap area during the study periods) to estimate the changes in vertical profiles, leaf area density (LAD), and canopy complexity.

Gap Detection and Classification
We set the thresholds of height and area for detecting open canopies as 5 m (which was almost half of the CHMs) and 10 m 2 (based on the gap area calculation from mobile LiDAR datasets and for convenience of calculation), respectively, and identified the open canopies (OCs) in each CHM. We defined gaps with canopy openings caused by fallen or dead trees as forest gaps (FGp) and those caused by civic usage as farmland, cemeteries, or trails as artificial gaps (AGp) (Figure 1).
To investigate the size distributions of FGp and AGp, we estimated the scaling exponent (λ) that describes the extent to which disturbances are clustered (i.e., lower λ means frequent large disturbances, while higher λ means frequent small disturbances) [49][50][51] using methods from Asner et al. [51].

Estimating Changes of Vertical Canopy Distribution and Canopy Complexity
The leaf area index (LAI, m 2 /m 2 ) and canopy complexity (Rumple index, m 2 /m 2 ) are well known to correlate with the productivity of forest ecosystems [52][53][54][55]. LAI is usually calculated as half the total leaf area per unit surface area [56,57]. We estimated the LAI and leaf area density (LAD, m 2 /m 3 ) using the "lidR" package in R software [58] according to the method of Bouvier et al. [52], who calculated the LAI based on the Beer-Lambert theory. The LAD was estimated by dividing the LAI values by dz (the profile of LAD). The canopy complexity was also estimated using the Rumple index function within the "lidR" package. Canopy complexity was calculated as a three-dimensional surface area divided by a two-dimensional surface area, which denotes the structural diversity of canopies [53][54][55]59]. Section 2.6 is summarized in Table 1 and Figure 4. Please see Figure 4 below

Gap Detection and Classification
We set the thresholds of height and area for detecting open canopies as 5 m (which was almost half of the CHMs) and 10 m 2 (based on the gap area calculation from mobile LiDAR datasets and for convenience of calculation), respectively, and identified the open canopies (OCs) in each CHM. We defined gaps with canopy openings caused by fallen or dead trees as forest gaps (FGp) and those caused by civic usage as farmland, cemeteries, or trails as artificial gaps (AGp) (Figure 1).
To investigate the size distributions of FGp and AGp, we estimated the scaling exponent (λ) that describes the extent to which disturbances are clustered (i.e., lower λ means frequent large disturbances, while higher λ means frequent small disturbances) [49][50][51] using methods from Asner et al. [51].

Estimating Changes of Vertical Canopy Distribution and Canopy Complexity
The leaf area index (LAI, m 2 /m 2 ) and canopy complexity (Rumple index, m 2 /m 2 ) are well known to correlate with the productivity of forest ecosystems [52][53][54][55]. LAI is usually calculated as half the total leaf area per unit surface area [56,57]. We estimated the LAI and leaf area density (LAD, m 2 /m 3 ) using the "lidR" package in R software [58] according to the method of Bouvier et al. [52], who calculated the LAI based on the Beer-Lambert theory. The LAD was estimated by dividing the LAI values by dz (the profile of LAD). The canopy complexity was also estimated using the Rumple index function within the "lidR" package. Canopy complexity was calculated as a three-dimensional surface area divided by a two-dimensional surface area, which denotes the structural diversity of canopies [53][54][55]59]. Section 2.6 is summarized in Table 1 and Figure 4.

Pixel and Hexagon Height Model-Based Change Detection
With change detection, we found that changes in the distribution of canopies were irregular year by year ( Figure 5), while the annual VGr rates remained relatively stable throughout the study period. There were no great growth rate differences according to either the time interval or the grid size ( Table 2). Figure 5a presents an example of the overall distributions of growth and damaged areas in 2012-2015. Figure 5b and 5c show the vertical changes over two years (i.e., one-year term) to explain the growth and damaged area distributions. Detected growth areas where the canopy height increased usually resulted from lateral growth rather than vertical growth in HCs, which could mean that taller canopy overlapping occurred frequently. In the case of low canopies, lateral growth areas were pervasive in annual change detection, while vertical growth areas were sustained over total periods. Because LC areas included bare lands, the height remained relatively unchanged.
Furthermore, canopy closure (LC to HC) tended to be driven by LGr rather than VGr, possibly because the rate of LGr among taller trees was higher than that of VGr among young trees at the study site (Figure 5d). Besides, we found that disturbance and canopy loss were widespread over the whole study period (Figure 5b and 5c).

Pixel and Hexagon Height Model-Based Change Detection
With change detection, we found that changes in the distribution of canopies were irregular year by year ( Figure 5), while the annual VGr rates remained relatively stable throughout the study period. There were no great growth rate differences according to either the time interval or the grid size (Table 2). Figure 5a presents an example of the overall distributions of growth and damaged areas in 2012-2015. Figure 5b,c show the vertical changes over two years (i.e., one-year term) to explain the growth and damaged area distributions. Detected growth areas where the canopy height increased usually resulted from lateral growth rather than vertical growth in HCs, which could mean that taller canopy overlapping occurred frequently. In the case of low canopies, lateral growth areas were pervasive in annual change detection, while vertical growth areas were sustained over total periods. Because LC areas included bare lands, the height remained relatively unchanged.
Furthermore, canopy closure (LC to HC) tended to be driven by LGr rather than VGr, possibly because the rate of LGr among taller trees was higher than that of VGr among young trees at the study site ( Figure 5d). Besides, we found that disturbance and canopy loss were widespread over the whole study period (Figure 5b,c).

Pixel and Hexagon Height Model-Based Change Detection
With change detection, we found that changes in the distribution of canopies were irregular year by year ( Figure 5), while the annual VGr rates remained relatively stable throughout the study period. There were no great growth rate differences according to either the time interval or the grid size ( Table 2). Figure 5a presents an example of the overall distributions of growth and damaged areas in 2012-2015. Figure 5b and 5c show the vertical changes over two years (i.e., one-year term) to explain the growth and damaged area distributions. Detected growth areas where the canopy height increased usually resulted from lateral growth rather than vertical growth in HCs, which could mean that taller canopy overlapping occurred frequently. In the case of low canopies, lateral growth areas were pervasive in annual change detection, while vertical growth areas were sustained over total periods. Because LC areas included bare lands, the height remained relatively unchanged.
Furthermore, canopy closure (LC to HC) tended to be driven by LGr rather than VGr, possibly because the rate of LGr among taller trees was higher than that of VGr among young trees at the study site (Figure 5d). Besides, we found that disturbance and canopy loss were widespread over the whole study period (Figure 5b and 5c).

Continuous One-Year Vertical Growth Area
Due to widespread canopy loss and the complexity of the dynamics in both canopies ( Figure 5), we tried to detect continuous growth areas. Table 3 shows the proportion of the area in which the canopy continuously grew upward at pixel (via CHMs) and plot (via hexagonal grids) levels. The proportion of overlapping growth areas with positive values (i.e., gradual vertical canopy growth) based on the height difference for 2012-2015 showed a 15.9% correspondence at the pixel scale and a 38.9% correspondence at the plot scale.

Continuous One-Year Vertical Growth Area
Due to widespread canopy loss and the complexity of the dynamics in both canopies ( Figure 5), we tried to detect continuous growth areas. Table 3 shows the proportion of the area in which the canopy continuously grew upward at pixel (via CHMs) and plot (via hexagonal grids) levels. The proportion of overlapping growth areas with positive values (i.e., gradual vertical canopy growth) based on the height difference for 2012-2015 showed a 15.9% correspondence at the pixel scale and a 38.9% correspondence at the plot scale.

Open Canopy Change Detection
For 2012-2015, the total number of FGps and the sum of the area decreased (Table 4). However, the mean FGp area each year did not decrease over time; this might indicate that some forest gaps continued to open and some were closed during the study periods. Moreover, total FGp areas (Sum in Table 4) abruptly decreased between 2013 and 2014, which was also present in Figure 5c, as growth distributions increased. Finally, even though the scaling exponent λ in 2013 decreased slightly, λ showed consistent values over the years. Forest gaps were classified into two subclasses according to whether they were closed in the study periods (e.g., [28]): existing gap (EGp) and closed gap (CGp) ( Table 5). The number of existing gaps (EGp) (i.e., FGp in the same location detected for all years) increased slightly from 374 to 391 due to gap coalescence [e.g., 28], while the mean area of EGp decreased by 19.11 m 2 (Table 5). Although the number of EGp patches increased slightly, the mean area and maximum area diminished because some gaps were divided into two or more gaps as a result of tree growth. Moreover, areas of 10.00-94.75 m 2 (mean: 16 ± 9 m 2 ) were closed for the whole study period (Table 5). Table 5. Existing gaps and closed gaps in forest gaps.

Forest Gap (FGp) Gap Continuation (Existing Gap, EGp) Closed Gap (CGp)
Year The difference in the number of existing gaps in each year might be due to gap coalescence (e.g., Vepakomma et al., [27]).
In the case of AGp, the mean area increased by about 10 m 2 while the number of patches decreased ( Table 6). The decrement in the number of patches could result from the LGr of the edge trees ( Table 5). The notable difference between the changes in forest gap and artificial openings is that artificial opening areas either maintained their position or expanded their area, while natural forest gaps showed more dynamics (i.e., faster canopy closure rates). Moreover, the scaling exponent λ showed lower values than that of forest gaps, which may indicate frequent larger openings in artificial gaps.

Changes in Vertical Canopy Structures in High Canopies and Open Canopies
The differences in changes to canopy structures by canopy class were more notable in the changes of vertical canopy distribution. In particular, FGp in OC showed both the largest LAI values each year and the largest differences in LAI (0.12 ± 0.62 m 2 /m 2 ) during the study periods, followed by HC and AGp, which could indicate the positive effects of FGp on foliage ( Figure 6). Furthermore, FGp also exhibited the greatest differences in canopy complexity for 2012-2015 (0.93 ± 4.15 m 2 /m 2 ) followed by HC and AGp (there was no significant difference between HCs and AGp), which possibly indicated that FGp had an important role in increasing canopy structural diversity (Figure 6b). Figure 7 shows the LAD of each year and differences in LAD for 2012-2015 (Figure 7 'class_dLAD') in FGp (subclassified as EGp and CGp), AGp, and HC at the study site. In the total study period (three-year interval), vertical foliage distribution decreased at a lower height and increased at a higher height for every canopy class. Particularly among high canopies, 12 m heights divided the LAD increment and decrement, while 5-10 m heights divided them in the low canopies. Furthermore, in the FGp canopy, the increment in LAD was greater than those in other canopy classes, while the differences in LAD among heights in the AGp did not show clear differences. However, looking at the annual changes in Figure 7, these changes were inconsistent throughout the study periods.

What Are the Differences between the Canopy Structural Changes Derived from Annual Change Detection and Three-Year Interval Change Detection?
We found that each annual vertical growth rate and total LAI showed almost consistent values, while distributions of changes in the canopies and leaf area profiles were irregular year by year.

What Are the Differences between the Canopy Structural Changes Derived from Annual Change Detection and Three-Year Interval Change Detection?
We found that each annual vertical growth rate and total LAI showed almost consistent values, while distributions of changes in the canopies and leaf area profiles were irregular year by year. Moreover, the distribution of the annual growth regions showed low correspondence with the three-year interval growth regions (Table 3). We speculated that this low correspondence and the irregular changes in the canopies might have resulted from frequent changes in urban canopies (i.e., dynamics), weather conditions during data acquisition, grid size, and the DOYs in data acquisition.

In Terms of the Aspects of the Urban Canopy Dynamics and the Day of the Year (DOY) of Data Acquisition
We found that the lateral canopy changes were irregular annually. By investigating Figure 5b and 5c, the first-year change detection in both HC and LC showed a higher distribution ratio for the damaged area than the others. In addition, the extent of the LGr areas was highest between 2013 and 2014, the second year (Figure 5b). Furthermore, Figure 6 shows that LAI and the canopy complexity (Rumple index) of all classes in 2013 were slightly lower than in previous years. These may show the importance of the DOYs of data acquisition for the one-year interval LiDAR survey. Since the DOY of data acquisition in 2013 was in the early growing season, while the DOY of data acquisition in 2014 was in the mid-growing season, the phenology might differ. This means that the point density of canopies in 2013 might have been relatively lower than that in other years, which could have resulted from the presence of more extensive damaged areas in the first year and smaller damaged areas in second year. However, considering the third year (the longest term of 401 d) when damaged areas occupied the largest area ( Figure 5) and changes in the leaf area density did not show clear differences (Figure 7), we could assume that both the DOY of data acquisition and other environmental variables might have affected the results.
Overall, the forest canopy appeared to change dynamically in each study year ( Figures 5 and 7). This may indicate that growth and disturbance events happened simultaneously and frequently in the study site. Table 7 shows the distribution of the yearly canopy changes in damaged area (DA) (HC to LC) and CaC (LC to HC). Row d in Table 7 shows the ratio of continuous change areas (open to closed or closed to open), which might indicate that dynamics of the canopy or uncertainty in the change detection using one-year interval LiDAR datasets.  Total (e: a + b + c − d) 100 100 Although the datasets for 2012-2015 were acquired on DOYs 283, 136, 167, and 303, respectively, the one-year interval comparisons showed the potential to estimate the VGr of the canopy surface, given the association between the annual VGr rate and DOY (Table 2) and to detect newly opened forest gaps during the study periods.

Regarding the Aspect of Grid Size
For the weather conditions, the canopy heights could have been influenced by windy conditions creating tree canopy movement. Furthermore, we were uncertain whether the trees had grown vertically. These could have resulted in low correspondences when detecting gradual growth areas at smaller grid sizes (Table 3), which could indicate that a smaller grid size could be easily affected by canopy movement or growth direction, while a larger grid size would not. In addition, as Roussel et al. [47] concluded, smaller plot sizes could cause larger errors in HMAX estimations; thus, the plot and grid size may have affected a higher degree of correspondence in the growth area for the plot-level survey.
However, in the case of detecting the gap opening areas, change detection with smaller grid sizes could be more useful than a larger size (Table 4, Figures 5 and 6). Since gap edge trees expand their foliage in the lateral direction ( Figures 5 and 7), we might not have seriously considered whether they have exactly grown.

What
Are the Characteristics of the Structural Changes According to the Different Canopy Classes (e.g., High Canopies and Low Canopies) in the Urban Forest?

Forest Gap Effects on Canopy Dynamics in Urban Forests
During the study period (2012)(2013)(2014)(2015), the overall distribution of the canopy growth area in HC was about 58.6% (VGr area: 17.1%, LGr area: 41.5%), while the damaged areas accounted for about 32.3% of the total area, which suggests active canopy growth and erosion in the study site. Meanwhile, about 26.2% of the LC area exhibited canopy growth, while the remaining 73.8% of the area was either damaged (14.6%) or unchanged (59.2%). The higher percentage of unchanged areas in LC probably resulted from areas of bare land or terrain not covered by canopies. Moreover, based on the VGr rate, both HCs and LCs seemed to grow upward at a steady rate (Table 2). However, the VGr rate of LCs was relatively lower than that of HCs; this difference might be due to civic usage and disturbances in LC.
AGp formed via anthropogenic activities (e.g., building graveyards, sports facilities, and recreation areas) has prolonged changes in LAI, canopy complexity, and LAD (Figures 6 and 7). The particularly slow or lower changes of AGp may have driven the low overall and VGr rates in LC. Regardless, the canopies of FGp showed the greatest changes in canopy complexity, LAI, and complexity ( Figures 6  and 7), exhibiting the well-known effects of FGp, which foster dynamic changes in the forest canopy [43]. Since canopy complexity and LAI are important parameters of forest health assessments, the FGp at the study site may have positively stimulated vegetation growth [55,56]. However, Figure 8 suggests that FGp could be also vulnerable to disturbances, agreeing with Vepakomma et al. [27].
In terms of the size distribution of physical disturbances or canopy openings, we could also find notable differences between FGp and AGp. In the case of FGp, the decrement of the scaling exponent in the first year (2012-2013) might result from frequent closures in smaller gaps rather than by opening larger areas (the maximum gap area and total area were similar between 2012-2013). From 2014 to 2015, it seemed that, due to forest gap closure, maximum and total gap area decreased, and the scaling exponent increased slightly. In the case of AGp, the decrement in the number of patches in the first year (2012-2013) and more frequent large openings might reflect slightly lower scaling exponent values in 2013 (the maximum gap area increased in 2013). From 2014 to 2015, it seemed that because of gap closure, maximum gap area was divided into several patches, and the scaling exponent increased slightly. Finally, since AGp had lower λ values, larger mean gap areas, and slower gap closure rates (6% closed in AGp, while 36% closed in FGp in all total periods, 2012-2015), we can assume that AGp was more vulnerable to physical disturbances; therefore, frequent larger openings occurred in AGp than in FGp [6,32].

Forest Gap Closure
The FGp exhibited a closure rate of 5.48 m 2 /year (Table 5), which appeared to be driven by LGr (Figure 5d) above a height of 5-10 m (i.e., HC) ( Figure 7). As shown in Figure 5d, LGr occupied most of the area in CaC (96.59%). Based on these results, the LGr rate of maturity consisted of the gap boundaries in trees that appeared to be greater than the VGr rate of the young trees or shorter vegetation on the floor. Many studies have reported that small disturbances such as branch and tree falls could be closed by the lateral extension of gap-edge trees [28,44,60,61]. At the present study site, FGps were smaller in size (mean area: 34-42 m 2 ) than AGp (mean area: 226-279 m 2 ). Considering a 5.48 m 2 /year rate of closure, FGp could close over short time scales. However, the larger gaps and reduced canopy dynamics in AGp areas led to more extensive gaps that remained open over long periods.

Conclusions
In this study, we estimated three-dimensional canopy changes in an urban forest using fouryear annual LiDAR datasets. Based on our results, we concluded that the vertical growth rates derived from one-year change detections showed similar values to three-year interval change detection. However, the regions in which the canopies grew up continuously differed according to the grid sizes (pixel-and plot-level grids), and they occupied only 15.9% and 38.9% of the growth areas detected in three-year interval comparisons. Moreover, the distributions in VGr, LGr, and damaged areas in the annual change detections showed irregular changes according to the year. These results might indicate that even though three-year interval (or larger) change detection could be a more efficient way to estimate the growth rates in a forest in terms of the costs and efforts, oneyear interval change detection with airborne LiDAR datasets would be useful for monitoring and understanding what areas are vulnerable to disturbance (three-year interval change detection might detect open canopies as closed canopies).
We also concluded that the leaf area index showed consistent values year by year, while the leaf area profile (e.g., leaf area density) showed inconsistent yearly changes in both high and low canopies. Among the low canopies, we could conclude that they were vulnerable to physical damage,

Forest Gap Closure
The FGp exhibited a closure rate of 5.48 m 2 /year (Table 5), which appeared to be driven by LGr (Figure 5d) above a height of 5-10 m (i.e., HC) ( Figure 7). As shown in Figure 5d, LGr occupied most of the area in CaC (96.59%). Based on these results, the LGr rate of maturity consisted of the gap boundaries in trees that appeared to be greater than the VGr rate of the young trees or shorter vegetation on the floor. Many studies have reported that small disturbances such as branch and tree falls could be closed by the lateral extension of gap-edge trees [28,44,60,61]. At the present study site, FGps were smaller in size (mean area: 34-42 m 2 ) than AGp (mean area: 226-279 m 2 ). Considering a 5.48 m 2 /year rate of closure, FGp could close over short time scales. However, the larger gaps and reduced canopy dynamics in AGp areas led to more extensive gaps that remained open over long periods.

Conclusions
In this study, we estimated three-dimensional canopy changes in an urban forest using four-year annual LiDAR datasets. Based on our results, we concluded that the vertical growth rates derived from one-year change detections showed similar values to three-year interval change detection. However, the regions in which the canopies grew up continuously differed according to the grid sizes (pixel-and plot-level grids), and they occupied only 15.9% and 38.9% of the growth areas detected in three-year interval comparisons. Moreover, the distributions in VGr, LGr, and damaged areas in the annual change detections showed irregular changes according to the year. These results might indicate that even though three-year interval (or larger) change detection could be a more efficient way to estimate the growth rates in a forest in terms of the costs and efforts, one-year interval change detection with airborne LiDAR datasets would be useful for monitoring and understanding what areas are vulnerable to disturbance (three-year interval change detection might detect open canopies as closed canopies).
We also concluded that the leaf area index showed consistent values year by year, while the leaf area profile (e.g., leaf area density) showed inconsistent yearly changes in both high and low canopies. Among the low canopies, we could conclude that they were vulnerable to physical damage, particularly in AGps. In the case of FGps, they also drove the growth and disturbance aspects of forest dynamics and served to diversify the canopy structure; most of them were closing in a lateral direction by gap-edge trees. Furthermore, high canopies expanded their foliage at a higher height (12 m) than those of low canopies (5 m).
Spatial-temporal distribution and gap characteristics could provide key information about what state these structures are in (e.g., close to natural forest state) vertically and horizontally through comparison with other forests. Moreover, based on the gap types (i.e., natural or human-induced openings), our results could guide management priorities, such as designing a conservation area or planning a restoration, by giving comprehensive information on which parts of a forest are vulnerable to disturbances or require better, more efficient management.
Although we found that one-year interval change detections showed the potential to achieve detailed information about canopy changes, further study is warranted. For example, since detecting changes in the distribution of canopies might be easily affected by the grid size and the DOY of data acquisition, an adequate grid size and DOY should be suggested. In particular, in the case of DOY, seasonal variations in the canopies over a year could be studied with a mobile LiDAR, and this may both hint at the DOY of data acquisition via airborne LiDAR survey and inform about seasonal changes in canopy structures.