Retrieval of Aerodynamic Parameters in Rubber Tree Forests Based on the Computer Simulation Technique and Terrestrial Laser Scanning Data

Rubber trees along the southeast coast of China always suffer severe damage from hurricanes. Quantitative assessments of the capacity for wind resistance of various rubber tree clones are currently lacking. We focus on a vulnerability assessment of rubber trees of different clones under wind disturbance impacts by employing multidisciplinary approaches incorporating scanned points, aerodynamics, machine learning and computer graphics. Point cloud data from two typical rubber trees belonging to different clones (PR107 and CATAS 7-20-59) were collected using terrestrial laser scanning, and a connection chain of tree skeletons was constructed using a clustering algorithm of machine learning. The concept of foliage clumps based on the trunk and first-order branches was first proposed to optimize rubber tree plot 3D modelling for simulating the wind field and assessing the wind-related parameters. The results from the obtained phenotypic traits show that the variable leaf area index and included angle between the branches and trunk result in variations in the topological structure and gap fraction of tree crowns, respectively, which are the major influencing factors relevant to the rubber tree’s capacity to resist hurricane strikes. The aerodynamics analysis showed that the maximum dynamic pressure, wind velocity and turbulent intensity of the wind-related parameters in rubber tree plots of clone PR107 (300 Pa, 30 m/s and 15%) are larger than that in rubber tree plots of clone CATAS-7-20-59 (120 Pa, 18 m/s and 5%), which results in a higher probability of local strong cyclone occurrence and a higher vulnerability to hurricane damage.


Introduction
Wind is a common disturbance agent in forest ecosystems. In recent years, the frequency and severity of wind damage through hurricanes have been increasing, causing significant losses in forestry [1]. Damage due to hurricanes has caused wind throw in forests and reduced the carbon As 3D laser scanners become more accurate and affordable and their ability to capture three-dimensional spatial arrangements increases, terrestrial laser scanning (TLS) is becoming an increasingly popular resource in forest management and ecology. Biological parameters extracted from TLS data, such as the gap fraction [27], tree volume [28], leaf area density [29] and tree diameter at breast height [30], can be successfully assessed without the need for arduous manual measurements. Furthermore, many algorithms for trees using TLS data have been derived, such as algorithms for individual leaf property extraction [31], leaf phenotypic characteristic calculation [32] and photosynthetic and non-photosynthetic part separation [33]. TLS is becoming increasingly useful for providing information to support forest observations, although the terms of the algorithm aspect for the detailed description of the whole tree from scanned data still require further improvement.
In summary, while many approaches can simulate wind damage in forests, models for hurricane damage on forests are still worthy of study. The following issues need to be addressed in the current researches. First, the forest stand is always considered to be porous media, and the internal morphological structure of the canopy, plant spacing or planting arrangement cannot be described accurately. Second, developing an individual tree scale model, including the phenotypic characteristics, the biological parameters and spatial information of vegetative elements in forest canopies, requires high computational complexity. Third, given the dangerous environment of a hurricane, using sensors to effectively measure the wind-related parameters (i.e., dynamic pressure, velocity and turbulent intensity) in various forest plots is infeasible due to heavy rain affection and unstable power supply.
We must therefore find an adequate compromise for modelling forest stands and propose the following: (1) Based on the point cloud data obtained by terrestrial laser scanning, a tree skeleton can be accurately identified and the trunk and multiple first-order branches automatically separated.
(2) A new concept of foliage clump composition based on trunks and different first-order branches is suggested, in which an optimized rubber tree canopy is modelled by several foliage clumps, and forest parameter retrievals of each foliage clump are performed to simplify the representation of the tree models and to depict the spatial structures of the trees as much as possible. (3) By applying aerodynamic modelling with 3D meshing for various forest plot scenarios composed of rubber tree models of different clones based on real planting conditions, here, we analyse the performance of rubber trees of different clones under a strong hurricane load.

Study Site and Data Collection
The study area was located in Hainan Dan Zhou (north-western area of Hainan Island, 109 • 430-109 • 510E; 19 • 280-19 • 380N) within a rubber tree plantation. The topography of the study area is a hilly plateau with an elevation of 188 m above sea level at the centre. The plateau is surrounded by flat lands with elevations of 20-160 m. As China's largest rubber production base, Hainan Island is continuously increasing the cultivation of rubber trees. The plantation has converted over 5000 ha of cultivated land and tropical rain forest since it was established in 1957. Of these lands, nearly 3000 ha were planted with natural rubber [34]. A wide variety of rubber trees clones was planted on the plantation, and in this study, two typical clones were selected: PR107 and CATAS7-20-59 ( Figure 1). Based on prior knowledge, CATAS7-20-59 is more resistant to hurricanes than is PR107, although a quantitative analysis has not been performed. Both rubber trees selected in this study belong to the Chinese Academy of Tropical Agricultural Sciences Proving Ground and are managed identically. In a previous study of rubber tree property retrieval [6], we proved that rubber trees belonging to the same clone have similar topological skeleton structures and phenotypic characteristics. Therefore, two typical rubber trees of each clone (PR107 and CATAS 7-20-59) were selected as the references for constructing the corresponding tree models. The point clouds of the two typical rubber trees, the canopy and skeleton, were obtained using a Leica C10 TLS system in May 2018. The Leica C10 instrument is a 532-nm phase-based scanner with a 360° × 270° upward field-of-view and a laser rate of 5000 points per second. The instrument setting height was 1.15 m, and the range measurement accuracy was ±1.5 mm at a distance of approximately 3.5 m. The circular laser beam diameter and beam divergence of the scanner exit were 3 mm and 0.22 mrad, respectively, yielding a minimum distance between consecutive beams of approximately 0.4 mm at a distance of 3.5 m from the instrument, with a scanning minimum deflection angle of 0.057°. Each tree plot was recorded using three scans around the target tree and using the "normal scanning" mode. The registration of the three point clouds corresponding to the three scans was carried out by manual adjustment based on a procedure in the software Cyclone (Leica Geosystem, Switzerland).
The structures of the two typical rubber trees represented by the TLS point clouds are illustrated in Figure 2a,b. The ground-based scanned point density values were 117,615 points (pts) m −2 and 88,249 pts m −2 for Rubber Tree 1 (PR107) and Rubber Tree 2 (CATAS7-20-59), respectively. The individual tree height, branch diameter, crown width, single-leaf area and angle between the branches were measured in situ. The tree height was measured using a Vertex IV hypsometer. The branch diameter was measured using a diameter tape. Crown widths were obtained as two values measured along two perpendicular directions from the treetop location. The single-leaf area and the angles between the trunk and branches were measured using an LI-3000C portable area metre and an angle measurement device, respectively. In a previous study of rubber tree property retrieval [6], we proved that rubber trees belonging to the same clone have similar topological skeleton structures and phenotypic characteristics. Therefore, two typical rubber trees of each clone (PR107 and CATAS 7-20-59) were selected as the references for constructing the corresponding tree models. The point clouds of the two typical rubber trees, the canopy and skeleton, were obtained using a Leica C10 TLS system in May 2018. The Leica C10 instrument is a 532-nm phase-based scanner with a 360 • × 270 • upward field-of-view and a laser rate of 5000 points per second. The instrument setting height was 1.15 m, and the range measurement accuracy was ±1.5 mm at a distance of approximately 3.5 m. The circular laser beam diameter and beam divergence of the scanner exit were 3 mm and 0.22 mrad, respectively, yielding a minimum distance between consecutive beams of approximately 0.4 mm at a distance of 3.5 m from the instrument, with a scanning minimum deflection angle of 0.057 • . Each tree plot was recorded using three scans around the target tree and using the "normal scanning" mode. The registration of the three point clouds corresponding to the three scans was carried out by manual adjustment based on a procedure in the software Cyclone (Leica Geosystem, Switzerland).
The structures of the two typical rubber trees represented by the TLS point clouds are illustrated in Figure 2a,b. The ground-based scanned point density values were 117,615 points (pts) m −2 and 88,249 pts m −2 for Rubber Tree 1 (PR107) and Rubber Tree 2 (CATAS7-20-59), respectively. The individual tree height, branch diameter, crown width, single-leaf area and angle between the branches were measured in situ. The tree height was measured using a Vertex IV hypsometer. The branch diameter was measured using a diameter tape. Crown widths were obtained as two values measured along two perpendicular directions from the treetop location. The single-leaf area and the angles between the trunk and branches were measured using an LI-3000C portable area metre and an angle measurement device, respectively.

Data Preprocessing
First, based on the wood-leaf separation method [33], LiDAR points were separated into branch point clouds (p ) and leaf point clouds (p ). The second step required filtering noise from the TLS point cloud, such as isolated laser points. The filtering noise algorithm [35] was adopted. We used balls with radii of 10 cm and removed all the balls containing less than 5 points. Figure 3 shows the wood-leaf separation results and noise reduction results. Separating the scanned points of the rubber trees into two parts (i.e., branch and leaf), where red represents branch point clouds ( ) and green represents leaf point clouds ( ). The segmentation results include a typical rubber tree of (a) clone PR107 and (b) clone CATAS 7-20-59.

Stratifying Branch Points and Obtaining the Central Points of Each Layer
According to the vertical heights of the individual trees, the branch point cloud is stratified into layers from the ground to the highest branch tip, and the number of layers depends on the tree height and the height interval ∆ℎ.
After obtaining the branch point clouds ( ) of every layer, we extract the central point from the branch point clouds ( ) in each layer using the cluster algorithm based on the Euclidean distance

Data Preprocessing
First, based on the wood-leaf separation method [33], LiDAR points were separated into branch point clouds (p b i ) and leaf point clouds (p l i ). The second step required filtering noise from the TLS point cloud, such as isolated laser points. The filtering noise algorithm [35] was adopted. We used balls with radii of 10 cm and removed all the balls containing less than 5 points. Figure 3 shows the wood-leaf separation results and noise reduction results.

Data Preprocessing
First, based on the wood-leaf separation method [33], LiDAR points were separated into branch point clouds (p ) and leaf point clouds (p ). The second step required filtering noise from the TLS point cloud, such as isolated laser points. The filtering noise algorithm [35] was adopted. We used balls with radii of 10 cm and removed all the balls containing less than 5 points. Figure 3 shows the wood-leaf separation results and noise reduction results. ) and green represents leaf point clouds ( ). The segmentation results include a typical rubber tree of (a) clone PR107 and (b) clone CATAS 7-20-59.

Stratifying Branch Points and Obtaining the Central Points of Each Layer
According to the vertical heights of the individual trees, the branch point cloud is stratified into layers from the ground to the highest branch tip, and the number of layers depends on the tree height and the height interval ∆ℎ.
After obtaining the branch point clouds ( ) of every layer, we extract the central point from the branch point clouds ( ) in each layer using the cluster algorithm based on the Euclidean distance

Stratifying Branch Points and Obtaining the Central Points of Each Layer
According to the vertical heights of the individual trees, the branch point cloud is stratified into LevelNum layers from the ground to the highest branch tip, and the number of layers depends on the tree height H t and the height interval ∆h.
Remote Sens. 2020, 12, 1318 6 of 22 After obtaining the branch point clouds (p b i ) of every layer, we extract the central point from the branch point clouds (p b i ) in each layer using the cluster algorithm based on the Euclidean distance [36] with a distance threshold of dist 2 . If the distance between each p b i in every layer is less than dist 1 , then the set of p b i is identified as one class in this case. Processed by our algorithm, the central point c l j c l j,x , c l j,y , c l j,z of every layer is extracted, where the jth central points in the lth layer are denoted by c l j . The central point connection between the adjacent layers is based on the minimum Dijkstra distance. As shown in Figure 4, three branch point layers and the central point of each layer are marked with different colours, and the black lines represent the connections between adjacent layers to produce the 3D tree skeleton.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 22 [36] with a distance threshold of . If the distance between each in every layer is less than , then the set of is identified as one class in this case. Processed by our algorithm, the central point , ， , ， , of every layer is extracted, where the j th central points in the th layer are denoted by . The central point connection between the adjacent layers is based on the minimum Dijkstra distance. As shown in Figure 4, three branch point layers and the central point of each layer are marked with different colours, and the black lines represent the connections between adjacent layers to produce the 3D tree skeleton.

Adopting the Cylinder Model to Form the Tree Skeleton
We define three categories of nodes in bottom-up order to highlight the central points that will be used. 1. Root node: root node is the lowermost central point.

Bifurcation node: bifurcation node
( ∈ ) has two or more connected child nodes. 3. Edge node: edge node ( ∈ ) has no connected child nodes (the nearest central point to is in the upper layer). In Figure 5, pink points represent the extracted root node , red points represent the extracted edge nodes and green points represent the extracted bifurcation nodes .

Adopting the Cylinder Model to Form the Tree Skeleton
We define three categories of nodes in bottom-up order to highlight the central points that will be used.

1.
Root node: root node c 1 1 is the lowermost central point.

2.
Bifurcation node: bifurcation node n l j (n l j ∈ c l j ) has two or more connected child nodes.

3.
Edge node: edge node e l i (e l i ∈ c l j ) has no connected child nodes (the nearest central point to c l j is in the upper layer).
In Figure 5, pink points represent the extracted root node c 1 1 , red points represent the extracted edge nodes e l j and green points represent the extracted bifurcation nodes n l j . Our algorithm begins from root node c 1 1 . According to the k-nearest neighbour and plant growth characteristics, the tree skeleton is transformed into the connection chains from the root node c 1 1 to every edge node e l i . The connection chain between each central point belonging to the adjacent layer is represented by a cylinder D l,l−1 i,j c l i , c l−1 j , r layer = 2, 3, 4, · · · · · · , levelnum , where r is the radius. According to the plant forest parameters, the diameter of the branch depends on the height; therefore, the lower the layer of the cylinder, the larger is the radius of the cylinder. For each cylinder D l,l−1 i,j c l i , c l−1 j , r , the radius is based on Dijkstra's route distance along the connection chain from the central point of the current cylinder to the root node, and this distance is denoted by Dis(c l i , c 1 1 ). The radius r of each cylinder can be obtained as follows: Remote Sens. 2020, 12, 1318 where λ is a constant related to the tree clones. As shown in Figure 6, the skeletons of Rubber Tree 1 (PR107) and Rubber Tree 2 (CATAS7-20-59) were reconstructed.

Adopting the Cylinder Model to Form the Tree Skeleton
We define three categories of nodes in bottom-up order to highlight the central points that will be used. 1. Root node: root node is the lowermost central point. 2. Bifurcation node: bifurcation node ( ∈ ) has two or more connected child nodes. 3. Edge node: edge node ( ∈ ) has no connected child nodes (the nearest central point to is in the upper layer). In Figure 5, pink points represent the extracted root node , red points represent the extracted edge nodes and green points represent the extracted bifurcation nodes . Our algorithm begins from root node . According to the k-nearest neighbour and plant growth characteristics, the tree skeleton is transformed into the connection chains from the root node to every edge node . The connection chain between each central point belonging to the adjacent layer is represented by a cylinder , , , , { = 2,3,4, ⋯ ⋯ , }, where is the radius. According to the plant forest parameters, the diameter of the branch depends on the height; therefore, the lower the layer of the cylinder, the larger is the radius of the cylinder. For each cylinder , , , , , the radius is based on Dijkstra's route distance along the connection chain from the central point of the current cylinder to the root node, and this distance is denoted by ( , ) . The radius r of each cylinder can be obtained as follows: where is a constant related to the tree clones. As shown in Figure 6, the skeletons of Rubber Tree 1 (PR107) and Rubber Tree 2 (CATAS7-20-59) were reconstructed.

Recognizing the Trunk and First-Order Branches
To more accurately describe the tree skeleton, we divided the tree skeleton into the trunk and many first-order branches. For the rubber tree, the trunk, at least near the base of the tree, is often nearly straight and extends with minimum changes in the growth angle. The trunk is a chain that must start from the root node (c 1 1 ), and the rest of the chains originating from the trunk will be defined as different first-order branches. Therefore, our algorithm is based on Dijkstra's algorithm and seeks the connection chain regarding the trunk depending on the rule of minimal change in the growth angle (see Appendix A). For the first-order branches, the connection composed of the central point is searched from the bifurcation nodes on the trunk to the other edge node. In Figure 7, we use different colours to represent our derivation results for the trunk and different first-order branches.

Recognizing the Trunk and First-Order Branches
To more accurately describe the tree skeleton, we divided the tree skeleton into the trunk and many first-order branches. For the rubber tree, the trunk, at least near the base of the tree, is often nearly straight and extends with minimum changes in the growth angle. The trunk is a chain that must start from the root node ( ), and the rest of the chains originating from the trunk will be defined as different first-order branches. Therefore, our algorithm is based on Dijkstra's algorithm and seeks the connection chain regarding the trunk depending on the rule of minimal change in the growth angle (see Appendix A). For the first-order branches, the connection composed of the central point is searched from the bifurcation nodes on the trunk to the other edge node. In Figure 7, we use different colours to represent our derivation results for the trunk and different first-order branches.

Determining Foliage Clumps based on the Trunk and First-Order Branches
According plant physiological theory, positive pressure resulting from metabolic pumping by the roots draws in water and nutrients from the soil to the leaves via sieve tubes and vessels, and similar hydraulic resistance is observed within vascular tissue via the same branch transport [37], which means that leaves belonging to the same first-order branch have not only a close spatial position and similar illumination but also almost equivalent nutrients and physiological properties [38]. We therefore defined a foliage clump as a leaf collection on the same trunk or first-order branch. The 3D watershed segmentation method [39] was employed to segment the foliage clumps based on the spatial locations of the trunk and first-order branches, and the segmented foliage clumps are represented by different colours (Figure 8).

Determining Foliage Clumps based on the Trunk and First-Order Branches
According plant physiological theory, positive pressure resulting from metabolic pumping by the roots draws in water and nutrients from the soil to the leaves via sieve tubes and vessels, and similar hydraulic resistance is observed within vascular tissue via the same branch transport [37], which means that leaves belonging to the same first-order branch have not only a close spatial position and similar illumination but also almost equivalent nutrients and physiological properties [38]. We therefore defined a foliage clump as a leaf collection on the same trunk or first-order branch. The 3D watershed segmentation method [39] was employed to segment the foliage clumps based on the spatial locations of the trunk and first-order branches, and the segmented foliage clumps are represented by different colours (Figure 8). Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 22

Retrieving the Foliage Clump Properties
After foliage clump extraction, the parameters of every foliage clump were calculated, including the leaf area index (LAI), crown volume, foliage clump volume, leaf area density and gap fraction. Thus, the 3D reconstruction of the total tree can be simplified into a model consisting of a tree branch skeleton and several foliage clumps. 1. LAI is the ratio of the total leaf area to ground area. First, after the single-leaf extraction and the number of leaves in each foliage clump were obtained using the method described in [29], Delaunay triangulation was adopted to deduce the area of each leaf. We acquired the LAI by computing the ratio of the sum of all leaf areas in each foliage clump to the projected area of each foliage clump. 2. Crown volume and foliage clump volume: A 3D convex hull algorithm [6] was used to deduce the tree crown volume and volume of each foliage clump. 3. Leaf area density: For each foliage clump, the leaf area density was expressed as the ratio of total leaf area to the volume of each foliage clump. 4. Gap fraction: The detailed derivation of the gap fraction of each foliage clump is available in Appendix B.

Retrieval of the Wind-Related Parameters in the Rubber Tree Plot
We constructed rubber tree models based on the foliage clump concept proposed in Section 2.3.4. Nine identical rubber tree models of each clone were placed into the plot according to the realistic line and row spacing distances (for PR107 and CATAS7-20-59, the line spacing and row spacing distances were 4 m and 9 m, respectively) to constitute the corresponding forest plot scenarios. After the 3D meshing for the two forest plot scenarios of each clone, the built forest scenario models were brought into the commercial computational fluid dynamics (CFD) software FLUENT (ANSYS Inc., Canonsburg, PA, USA) to simulate a full two-way interaction between the hurricane and forest stands and to calculate the wind-related parameters in the tree plots. The computational model using FLUENT software is based on a standard k-ϵ two equation model [40] (Appendix C affording the full description). The rubber tree plantation in Hainan usually suffers from Level 7 hurricanes with wind speeds of 13.9-17.1 m/s; therefore, we set the wind velocity to 15 m/s, and the results of the specific wind-related parameters in the forest plot, such as the velocity, dynamic pressure and turbulent intensity, are shown in the Results section.

Retrieving the Foliage Clump Properties
After foliage clump extraction, the parameters of every foliage clump were calculated, including the leaf area index (LAI), crown volume, foliage clump volume, leaf area density and gap fraction. Thus, the 3D reconstruction of the total tree can be simplified into a model consisting of a tree branch skeleton and several foliage clumps.

1.
LAI is the ratio of the total leaf area to ground area. First, after the single-leaf extraction and the number of leaves in each foliage clump were obtained using the method described in [29], Delaunay triangulation was adopted to deduce the area of each leaf. We acquired the LAI by computing the ratio of the sum of all leaf areas in each foliage clump to the projected area of each foliage clump.

2.
Crown volume and foliage clump volume: A 3D convex hull algorithm [6] was used to deduce the tree crown volume and volume of each foliage clump.

3.
Leaf area density: For each foliage clump, the leaf area density was expressed as the ratio of total leaf area to the volume of each foliage clump. 4.
Gap fraction: The detailed derivation of the gap fraction of each foliage clump is available in Appendix B.

Retrieval of the Wind-Related Parameters in the Rubber Tree Plot
We constructed rubber tree models based on the foliage clump concept proposed in Section 2.3.4. Nine identical rubber tree models of each clone were placed into the plot according to the realistic line and row spacing distances (for PR107 and CATAS7-20-59, the line spacing and row spacing distances were 4 m and 9 m, respectively) to constitute the corresponding forest plot scenarios. After the 3D meshing for the two forest plot scenarios of each clone, the built forest scenario models were brought into the commercial computational fluid dynamics (CFD) software FLUENT (ANSYS Inc., Canonsburg, PA, USA) to simulate a full two-way interaction between the hurricane and forest stands and to calculate the wind-related parameters in the tree plots. The computational model using FLUENT software is based on a standard k-two equation model [40] (Appendix C affording the full description). The rubber tree plantation in Hainan usually suffers from Level 7 hurricanes with wind speeds of 13.9-17.1 m/s; therefore, we set the wind velocity to 15 m/s, and the results of the specific wind-related parameters in the forest plot, such as the velocity, dynamic pressure and turbulent intensity, are shown in the Results section.

Properties of the Two Rubber Trees
In our algorithm, for Tree 1 (PR107), the ∆h (in Section 2.3.1.) and λ (in Section 2.3.2.) values were set as 0.05 m and 0.0778, respectively. For Tree 2 (CATAS7-20-59), the ∆h and λ values were set as 0.03 m and 0.0714, respectively.
In this study, some of the forest parameters obtained directly from the point cloud data are listed in Table 1. Other essential parameters were validated by in situ measurements and are listed in Table 2, where A, B, D and E represent first-order branches and C represents the trunk, as shown in Figure 7.   Note: E-W: east-west direction; N-S: north-south direction; C: trunk.
The growth properties of five foliage clumps of Tree 1 (PR107) and Tree 2 (CATAS7-20-59) are shown in Figure 8. We computed the total number of point clouds on each foliage clump, foliage clump volume, projection area, total number of leaves, leaf area, LAI, leaf area density and gap fraction and compared the parameters for the two rubber trees (Table 3).

Reconstruction of the Forest Plot Model
Nine typical and identical tree models belonging to each clone (PR107 or CATAS 7-20-59) were reconstructed according to the derived foliage clumps and branch attributes using our method, which constitute the corresponding forest plots of two clones. Then, the meshing step for the two forest plots was performed to facilitate the discrete solutions of the aerodynamic formulas (Equations (A9), (A12) and (A13) in Appendix C) for calculation of the aerodynamic parameters in the forest plots, as shown in Figure 9.

Reconstruction of the Forest Plot Model
Nine typical and identical tree models belonging to each clone (PR107 or CATAS 7-20-59) were reconstructed according to the derived foliage clumps and branch attributes using our method, which constitute the corresponding forest plots of two clones. Then, the meshing step for the two forest plots was performed to facilitate the discrete solutions of the aerodynamic formulas (Equations (A9), (A12) and (A13) in Appendix C) for calculation of the aerodynamic parameters in the forest plots, as shown in Figure 9.    11~12 m). In contrast, inside Forest Plot 2, the velocity is relatively evenly distributed and less than that of Forest Plot 1. Large speed changes occur on both sides of the rubber trees in rows (Y = 10~12 m and 19~22 m) because the hurricane propagation meets the barrier composed of rubber trees and rushes through the spacing between the trees. Figure 11a,b presents the vertical velocity profiles (Y = 15 m). The maximal values exceed 24 m/s and are mainly located in the windward area above the tree crown of Forest Plot 1.  First, by comparing the horizontal section (Z = 7 m) of the velocity distribution for the two forest plots, Figure 10a,b shows that the velocity inside Forest Plot 1 has a larger fluctuation range and a speed of up to 24 m/s in the canopies of the second and third trees (approximately X = 7~8 m and 11~12 m). In contrast, inside Forest Plot 2, the velocity is relatively evenly distributed and less than that of Forest Plot 1. Large speed changes occur on both sides of the rubber trees in rows (Y = 10~12 m and 19~22 m) because the hurricane propagation meets the barrier composed of rubber trees and rushes through the spacing between the trees. Figure 11a,b presents the vertical velocity profiles (Y = 15 m). The maximal values exceed 24 m/s and are mainly located in the windward area above the tree crown of Forest Plot 1.

Analysing the Wind-Related Parameters in Forest Plots of Different Clones
Second, from Figure 10c,d, for Forest Plot 1 (PR107), the dynamic pressure (i.e., the kinetic energy per unit volume of a fluid particle) value is relatively higher than that of Forest Plot 2 (CATAS7-20-59) and changes drastically inside a single tree, especially in the canopy of the third tree (X = 11~13 m), reaching above 180 Pa. In contrast, in Forest Plot 2, the dynamic pressure tends to be stable and generally less than 120 Pa. Then, we focused on a vertical section of dynamic pressure (Figure 11c,d). A great change in pressure occurs near the windward area of the forest tops in Forest Plot 1 (Z = 16~42 m), and this change is not as marked in Forest Plot 2. Local maxima in the horizontal and vertical profiles of dynamic pressure are both formed in Forest Plot 1 (Figures 10c and 11c), with one on the leeside of the tree and another above the windward side of the crown, and the dynamic pressure in both places exceeds 180 Pa. Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 22 Second, from Figure 10c,d, for Forest Plot 1 (PR107), the dynamic pressure (i.e., the kinetic energy per unit volume of a fluid particle) value is relatively higher than that of Forest Plot 2 (CATAS7-20-59) and changes drastically inside a single tree, especially in the canopy of the third tree (X = 11~13 m), reaching above 180 Pa. In contrast, in Forest Plot 2, the dynamic pressure tends to be stable and generally less than 120 Pa. Then, we focused on a vertical section of dynamic pressure (Figure 11c,d). A great change in pressure occurs near the windward area of the forest tops in Forest Plot 1 (Z = 16~42 m), and this change is not as marked in Forest Plot 2. Local maxima in the horizontal and vertical profiles of dynamic pressure are both formed in Forest Plot 1 (Figures 10c and 11c), with one on the leeside of the tree and another above the windward side of the crown, and the dynamic pressure in both places exceeds 180 Pa.
Finally, Figure 10e,f shows the turbulent intensity (i.e., the ratio of the standard deviation of the fluctuating wind velocity to the mean wind speed) between the two forest plots. Turbulent intensity is a measure of the degree of pulsation of airflow velocity, and generally divided into three categories according to magnitude: high turbulence (5%-20%), medium turbulence (1%-5%) and low turbulence (less than 1%). The larger the turbulent intensity, the more chaotic the behaviour of the airflow motion. Upon combination with Figure 11e, it is obvious that unstable turbulence (up to 15%) occurs in both sides of Forest Plot 1, which means that the air flows in Forest Plot 1 present a high variation in wind velocity and generate a high wind shear force. Coupled with high LAI and larger crown volumes leading to the increment in the frontal area opposing wind flow and small gaps among the forests, hurricanes in Forest Plot 1 generate strong interactions with vegetative elements and are Finally, Figure 10e,f shows the turbulent intensity (i.e., the ratio of the standard deviation of the fluctuating wind velocity to the mean wind speed) between the two forest plots. Turbulent intensity is a measure of the degree of pulsation of airflow velocity, and generally divided into three categories according to magnitude: high turbulence (5%-20%), medium turbulence (1%-5%) and low turbulence (less than 1%). The larger the turbulent intensity, the more chaotic the behaviour of the airflow motion. Upon combination with Figure 11e, it is obvious that unstable turbulence (up to 15%) occurs in both sides of Forest Plot 1, which means that the air flows in Forest Plot 1 present a high variation in wind velocity and generate a high wind shear force. Coupled with high LAI and larger crown volumes leading to the increment in the frontal area opposing wind flow and small gaps among the forests, hurricanes in Forest Plot 1 generate strong interactions with vegetative elements and are prone to generate branch breakage and serious defoliation. However, in Forest Plot 2 (Figures 10f and 11f), the turbulent intensity changes are not obvious and are generally lower than 5%, hinting at minor fluctuations of wind velocity existing and air flowing more smoothly. Overall, the influence of inhomogeneities in the aerodynamic parameters of Forest Plot 1 is more pronounced than that of Forest Plot 2, which leads to a greater possibility of wind-induced tree body damage in Forest Plot 1.
We extracted a cuboid region (length x = 0~24 m, width y = 13~18 m and height z = 8~16 m in Forest Plots 1 and 2) to represent the specific wind-related parameter distribution in the middle of the forest canopy, as shown in Figure 12a,b. The range area graphs of velocity, dynamic pressure and turbulent intensity in the cuboids of the two forest plots are shown in Figure 12c,e. Seen from the graph of velocity (Figure 12c), the wind velocity in Plot 1 and Plot 2 tend to decrease when meeting the tree canopy due to the interception of the vegetative elements in the tree canopy. Meanwhile, Figure 12c shows a larger velocity fluctuation in the tree canopies at X = 6~8 m and 10~12 m of Forest Plot 1 than that in the Forest Plot 2. Figure 12d shows the distribution of dynamic pressure in the cuboid regions of the forest plots, and it can be generally seen that the trend of the dynamic pressure distribution is the same as that of the velocity distribution. The explanation of this case is the square of the velocity proportional to the dynamic pressure based on Bernoulli's equation [41]. Figure 12e shows the distribution of turbulent intensity in each cuboid region of the two forest plots. The turbulent intensity represents the intensity of wind velocity fluctuation [42]. A low gap fraction of the typical trees in Forest Plot 1 induced by higher leaf area density results in a local strong wind current entering the tree canopy, which causes the velocity to fluctuate greatly in the red area of the forest canopy at X = 3~4 m, 7~8 m and 12~13 m. Overall, the velocity, dynamic pressure and turbulent intensity of Forest Plot 1 are generally higher than those of Forest Plot 2, resulting in higher vulnerability to hurricane damage for Forest Plot 1 than for Forest Plot 2.
inhomogeneities in the aerodynamic parameters of Forest Plot 1 is more pronounced than that of Forest Plot 2, which leads to a greater possibility of wind-induced tree body damage in Forest Plot 1.
We extracted a cuboid region (length x = 0~24 m, width y = 13~18 m and height z = 8~16 m in Forest Plots 1 and 2) to represent the specific wind-related parameter distribution in the middle of the forest canopy, as shown in Figure 12a,b. The range area graphs of velocity, dynamic pressure and turbulent intensity in the cuboids of the two forest plots are shown in Figure 12c,e. Seen from the graph of velocity (Figure 12c), the wind velocity in Plot 1 and Plot 2 tend to decrease when meeting the tree canopy due to the interception of the vegetative elements in the tree canopy. Meanwhile, Figure 12c shows a larger velocity fluctuation in the tree canopies at X = 6~8 m and 10~12 m of Forest Plot 1 than that in the Forest Plot 2. Figure 12d shows the distribution of dynamic pressure in the cuboid regions of the forest plots, and it can be generally seen that the trend of the dynamic pressure distribution is the same as that of the velocity distribution. The explanation of this case is the square of the velocity proportional to the dynamic pressure based on Bernoulli's equation [41]. Figure 12e shows the distribution of turbulent intensity in each cuboid region of the two forest plots. The turbulent intensity represents the intensity of wind velocity fluctuation [42]. A low gap fraction of the typical trees in Forest Plot 1 induced by higher leaf area density results in a local strong wind current entering the tree canopy, which causes the velocity to fluctuate greatly in the red area of the forest canopy at X = 3~4 m, 7~8 m and 12~13 m. Overall, the velocity, dynamic pressure and turbulent intensity of Forest Plot 1 are generally higher than those of Forest Plot 2, resulting in higher vulnerability to hurricane damage for Forest Plot 1 than for Forest Plot 2.

Discussion
In this study, we adopted a balance between the forest scale and individual tree scale, used the concept of foliage clumps to facilitate rubber tree modelling and retained the original morphological characteristics of rubber trees. The simplified model developed in this paper could represent the first step towards a new generation of medium-view models for simulating broad-leaf trees. Simplification into foliage clumps decreases the number of simulated leaves for which calculations must be performed and retains the local phenotypic characteristics of the rubber tree crown. Multidisciplinary approaches from aerodynamics [43], forestry and computer graphics were combined to reconstruct various forest plot scenarios and to simulate the wind fields in different rubber tree plots to quantitatively analyse the distribution of the wind-related parameters. The method used in this paper is flexible and can be extended to any broad-leaved tree species because the branching habits of broad-leaved trees and their phenotypic characteristics are similar. Through computer simulation techniques, spatial planning, wind speeds and the phenotypic characteristics of tree species (such as the gap fraction, LAI, canopy volume, etc.) can be adjusted to quantitatively investigate the

Discussion
In this study, we adopted a balance between the forest scale and individual tree scale, used the concept of foliage clumps to facilitate rubber tree modelling and retained the original morphological characteristics of rubber trees. The simplified model developed in this paper could represent the first step towards a new generation of medium-view models for simulating broad-leaf trees. Simplification into foliage clumps decreases the number of simulated leaves for which calculations must be performed and retains the local phenotypic characteristics of the rubber tree crown. Multi-disciplinary approaches from aerodynamics [43], forestry and computer graphics were combined to reconstruct various forest plot scenarios and to simulate the wind fields in different rubber tree plots to quantitatively analyse the distribution of the wind-related parameters. The method used in this paper is flexible and can be extended to any broad-leaved tree species because the branching habits of broad-leaved trees and their phenotypic characteristics are similar. Through computer simulation techniques, spatial planning, wind speeds and the phenotypic characteristics of tree species (such as the gap fraction, LAI, canopy volume, etc.) can be adjusted to quantitatively investigate the aerodynamic parameters of different stands under windy conditions. Moreover, comparing the damage risk of rubber tree plots in this research provides a better understanding of the differences in damage patterns among the two rubber tree clones and provides guidelines for suitable silvicultural management practices.
The simulation results are difficult to verify empirically. Experimental manipulations of wind at the scale of forest stands are impossible. Obtaining measurements during a hurricane is dangerous, and it is difficult to protect sensitive equipment from wind and rain damage during extreme weather. Our approach to stand reconstruction is generally applicable to broad-leaved trees but will not transfer as readily to conifers, which remain poorly covered by existing wind risk models. This is due to several limitations. Terrestrial laser scanning is an inappropriate tool with which to detect conifer needles, as they are generally smaller than the characteristic laser spot size. It is difficult to represent needles using point clouds, and the relationship between needles and the mechanical load of wind is hard to predict.
Our simulation results reveal a contrast between Forest Plot 1 (PR107) and Forest Plot 2 (CATAS7-20-59). Trees in Forest Plot 1 become more vulnerable to damage than do those in Forest Plot 2 because of the stronger velocity, dynamic pressure and turbulent intensity, particularly between the trunk and first-order branches. The analytical conclusion from our algorithm is consistent with the experience of rubber tree growers. Tree characteristics, such as gap fraction and shape, may be significant predictors of wind damage. Rubber trees of clone PR107 present a higher LAI and lower gap fraction than do rubber trees of clone CATAS 7-20-59, which means that the leaves of the rubber trees of clone PR107 are denser and cause considerable wind drag on the tree crown. Previous studies of wind damage in forests have suggested that trees with large, dense crowns produce large wind drag, resulting in tree lodging [44,45]. Moreover, the rubber trees of clone CATAS 7-20-59, which present a small angle between the trunk and the first-order branches to form an ascending vase shape for the tree crown, have a considerably more stable tree crown structure than do the rubber trees of clone PR107, with a spread-out tree crown under the impact of wind. The effect of tree shapes on wind firmness is also well documented [44,45], and tree architecture itself can alter the magnitude and the spatial distribution of wind loading [46].

Conclusions
Based on our findings obtained from the simulation program, the following strategies are suggested to reduce the risk of windthrow:

1.
Trees with large or dense crowns are more vulnerable to windthrow than are trees with smaller open crowns. Crown modification techniques, such as pruning and topping to reduce the effective crown size and density, can considerably reduce the risk of windthrow. Where possible, creating gaps that are too large and exposing individual branches or foliage clumps through these types of cuts should be avoided.

2.
A wide variety of rubber tree clones is planted in the coastal areas of China. The choice of rubber tree clones should take into account the probability of wind damage. Before extensively promoting a new clone of rubber trees, our method can be used to analyse the forest parameters, determine their aerodynamic parameters under windy conditions and measure the resistance capability of tree clones.

3.
Quantification of wind damage under different forest cultivation practices (e.g., adjusting the spatial distance between trees or changing the arrangement of trees) in the forest can be explored using our method to analyse to identify suitable silvicultural management strategies for different rubber tree clones.
Such management strategies will promote wind-resistant tree characteristics and help mitigate the risk of tree failures and subsequent economic damage.

Acknowledgments:
The authors gratefully acknowledge the foresters in the Chinese Academy of Tropical Agricultural Sciences Proving Ground for their assistance with data collection and sharing their rich knowledge and working experience.

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

Appendix A. Rule of Minimal Change in the Growth Angle
To determine the trunk chain, first, beginning from the root node (c 1 1 ), the node connected to the root node (c 1 1 ) is stored in a queue until meeting a bifurcation node (n layer j ). Secondly, as suggested by Figure A1 constitutes trunk k = 1, 2, 3 · · · ; i = 1, 2, 3 · · · (A1) Finally, when the edge node is determined, the central point belonging to the chain that begins from the root node is classified as the trunk.

Acknowledgements:
The authors gratefully acknowledge the foresters in the Chinese Academy of Tropical Agricultural Sciences Proving Ground for their assistance with data collection and sharing their rich knowledge and working experience.

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

Appendix A. Rule of Minimal Change in the Growth Angle.
To determine the trunk chain, first, beginning from the root node ( ), the node connected to the root node ( ) is stored in a queue until meeting a bifurcation node ( ). Secondly, as suggested by Figure A1, to implement the trunk chain determination, our algorithm presents a minimal change judgement method based on an angle value comparison. For bifurcation nodes ( ) on the trunk chain with multiple central points ( and ) in the upper layer, we define the included angle ∠ ， ⃗ ， ， ⃗ as . Similarly, we define the second included angle ∠ ， ⃗ ， ， ⃗ as . Thirdly, the chain that belongs to the trunk is determined by comparing and . We choose the chain with the minimum change in the direction angle as the trunk chain, and the angle comparison equation is as follows: Finally, when the edge node is determined, the central point belonging to the chain that begins from the root node is classified as the trunk.

Appendix B. Derivation of the Gap Fraction.
Each foliage clump is a collection of leaves and records the partial features of the whole tree crown. According to theoretical principles [47,48], by assuming that the inside of each foliage clump is a turbid medium, we have the following equation:

Appendix C. Standard k-Two Equation Model
Appendix C. 1

. Momentum Model
According to the momentum conservation law and Navier-Stokes equations, the momentum model is described by Equation (A9).
Here, ρ = 1.293 kg/m 3 and represents the density of a fluid; µ t (in kg/m) is the dynamic viscosity described by Equation (A16); p is the pressure of the wind; i and j represent any two directions of x, y and z in a three-dimensional coordinate system; x i and x j are the coordinate components in the i-direction and j-direction, respectively; u i (in m/s) is the wind speed in the i-direction; u j (in m/s) is the wind speed in the j-direction; and S i is the source term in the i-direction described by Equation (A10): where the LAI can be obtained in Section 2.3.5.

Appendix C.2. k − Model
The k − model overcomes the problem with flow close to the surface elements due to a singularity in the governing equations, and the k − model is described using the turbulent kinetic energy k and dissipation rate by Equations (A12) and (A13), respectively. The turbulent kinetic energy k and the dissipation rate in a fully developed neutral atmospheric boundary layer are defined by Equations (A14) and (A15), respectively.
As shown in Table A1, C µ , C ε1 , C ε2 , σ k and σ ε are momentum and turbulence equation constants used for simulating the forest canopy to retrieve some wind parameters. By combining Equations (A9), (A12) and (A13), the velocity, dynamic pressure and turbulent intensity can be deduced.