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

Abstract

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.

1. 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 storage potential of forests [2], possibly resulting in climatic changes [3]. Rubber trees (Hevea brasiliensis Müll.Arg.) are among the most important commercial sources of high-quality natural rubber in the world [4], which makes rubber trees one of the most important industrial crops for tropical regions in South and Southeast Asia [5]. In Hainan, China, a large number of rubber tree clones have been planted, although the rubber forests are subjected to periodic hurricanes. For example, in 2011 and 2014, the losses of tapping rubber trees reached 3.5 million and 4.7 million Yuan on the Leizhou Peninsula due to hurricanes Nasha and Rammasun, respectively. In 2016, tropical storm Dianmu (August 18th), severe tropical storm Mirinae (July 26th) and super hurricane Sarika (October 18th) caused extensive wind damage in the rubber plantations and led to relatively low temperatures [6]. To better predict wind damage to rubber trees, it is necessary to qualitatively and quantitatively measure the wind resistance capability of rubber trees of different clones under high-wind conditions [7].

Previous analyses have shown that one of the key issues in predicting and mitigating hurricane damage to trees is understanding the reasons behind tree damage [8,9]. Currently, the main methods for studying hurricane-induced tree damage are performed either at the large scale, forest stand scale or individual tree scale.

For the large scale, remote sensing data, such as satellite imagery, have been combined with meteorological data to assess and analyse forest damage before and after a hurricane, thereby providing guidance for forest management [10]. Commonly used satellites include the Google Earth satellite [11], Ice Cloud and Elevation Satellite (ICE Sat) [12], Landsat Thematic Mapper (TM) and SPOT 16-m [13], and these data have been applied to study wind-induced damage to forests in the Philippines, the United States and Finland, respectively. Furthermore, high-spatial resolution (1-m panchromatic and 4-m multispectral) IKONOS satellite imagery, which is finer in comparison to Landsat or SPOT 16-m imagery, has been used to reflect the disturbance severity of windstorm damage in north-eastern Minnesota [14]. These approaches that use satellite imagery have the advantage of being simple and therefore easily applicable; however, these approaches are limited in their power to retrieve detailed and specific wind-related parameters in forest canopies.

At the forest stand scale, the common approach includes linking risk models with simulation models of forest stand development, which has enabled the assessment of interactions between wind and trees (e.g., through wind-induced changes in stand structure, which in turn influence future susceptibility to disturbances) and the impacts to specific ecosystem services [15]. Simulation models of forest stands are divided into simple reconstructions [16] and accurate reconstructions [17]. Simple reconstructions simulate the canopy with a uniform and fixed shape for different tree species and improves upon the lack of preciseness in the tree model [18]. For example, Forest GALES effectively models the dynamics of wind damage to forest stands by treating the whole forest stand as a poroelastic continuous medium [19]. However, these models are currently limited to predictions for structurally uniform single-species stands [20]. Field measurements for modelling forest stands yield accurate reconstructions but are extremely time consuming and computationally expensive [21].

A few finer individual tree models have been developed. The two major research objects are the branch and the leaf. For the branch, measurements taken from trees growing in exposed and sheltered areas within two structurally similar forests have been used to investigate the influence of wind on the branch [22]. A multimodal approach was developed to represent the dynamic parameters of branches on trees during wind using complex models and finite element analyses [23]. Chiba [24] studied a hurricane-damaged sugi (Cryptomeria japonica D. Don) forest plantation by simulating physical stress changes in the tree stems caused by wind. Others have investigated the effects of wind on leaves. Models of tree geometry and leaf deformation have been combined to explore the role of wind mechanical parameters on leaf inclination angle distributions [25]. Zhu and Shao et al. [26] tested steady tulip tree (Liriodendron tulipifera L.) leaves under suspended conditions and obtained five critical wind speeds that can lead to leaf vibration. However, these models are either unable to simulate the wind parameters of the whole tree or can target only one type of tree.

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.

2. Materials and Methods

2.1. 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⁻² and 88,249 pts m⁻² 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.

2.2. Data Preprocessing

First, based on the wood-leaf separation method [33], LiDAR points were separated into branch point clouds (${\mathrm{p}}_{\mathrm{i}}^{\mathrm{b}}$) and leaf point clouds (${\mathrm{p}}_{\mathrm{i}}^{\mathrm{l}}$). 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.

2.3. Branch Skeleton Reconstruction

2.3.1. 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 $\Delta h.$

$$LevelNum=\frac{H_t}{\Delta h}$$

(1)

After obtaining the branch point clouds ($p_i^b$) of every layer, we extract the central point from the branch point clouds ($p_i^b$) in each layer using the cluster algorithm based on the Euclidean distance [36] with a distance threshold of $dis{t}_2$. If the distance between each $p_i^b$ in every layer is less than $dis{t}_1$, then the set of $p_i^b$ is identified as one class in this case. Processed by our algorithm, the central point $c_j^l\left(c_{j,x}^l,c_{j,y}^l,c_{j,z}^l\right)$ of every layer is extracted, where the jth central points in the lth layer are denoted by $c_j^l$. 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.

2.3.2. 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.

• Root node: root node $c_1^1$ is the lowermost central point.
• Bifurcation node: bifurcation node $n_j^l$ ($n_j^l\in c_j^l$) has two or more connected child nodes.
• Edge node: edge node $e_i^l(e_i^l\in c_j^l)$ has no connected child nodes (the nearest central point to $c_j^l$ 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_j^l$ and green points represent the extracted bifurcation nodes $n_j^l$.

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_i^l$. The connection chain between each central point belonging to the adjacent layer is represented by a cylinder $D_{i,j}^{l,l-1}\left(c_i^l,c_j^{l-1},r\right)\left\{layer=2,3,4,\cdots \cdots ,levelnum\right\}$, 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_{i,j}^{l,l-1}\left(c_i^l,c_j^{l-1},r\right)$, 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_i^l,c_1^1)$. The radius r of each cylinder can be obtained as follows:

$$r_{D_{i,j}^{l,l-1}}=\lambda \times \sqrt[3]{Dis\left(c_i^l,c_1^1\right)}$$

(2)

where $\lambda$ 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.

2.3.3. 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.

2.3.4. 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).

2.3.5. 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.

• 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.
• 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.
• 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.
• Gap fraction: The detailed derivation of the gap fraction of each foliage clump is available in Appendix B.

2.3.6. 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.

3. Results

3.1. 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. Studied rubber tree parameters, which were derived directly from the data.

Tree	Total Number of Leaf Points	Tree Crown Volume (m3)	Average Single-Leaf Area (cm2)	Tree Crown Projection Area (m2)	LAI (m2/m2)
Tree 1	117,615	168.73	75.53	16.08	3.08
Tree 2	88,249	142.51	79.54	12.65	2.62

. Other essential parameters were validated by in situ measurements and are listed in

Table 2. Forest parameters obtained from our method versus field measurements.

Tree	Height (m)/ Crown Diameter (m) (E-W) × (N-S)	Branch Diameter (cm) (Our Method/ Field Measurement)	Angle between the Trunk and the First-Order Branch (°) (Our Method/Field Measurement)
			$\angle (\mathbf{A},\mathbf{C})$	$\angle (\mathbf{B},\mathbf{C})$	$\angle (\mathbf{D},\mathbf{C})$	$\angle (\mathbf{E},\mathbf{C})$
Tree 1	15.36/ 3.85 × 5.71	A:21.6/22.1 B:22.3/20.8 C:28.7/30.5 D:25.3/23.9 E:18.7/20.8	45.19/ 47.23	53.14/ 49.36	47.37/ 45.64	60.72/ 57.56
Tree 2	17.13/ 3.07 × 5.59	A:20.7/22.8 B:16.4/15.7 C:35.1/36.8 D:25.6/27.3 E:18.6/19.5	42.36/ 41.78	37.89/ 40.25	34.47/ 32.92	43.91/ 42.24

Note: E-W: east–west direction; N-S: north–south direction; C: trunk.

, where A, B, D and E represent first-order branches and C represents the trunk, as shown in Figure 7.

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. The calculated growth properties of each foliage clump belonging to various branches.

Tree	Foliage Clump Belonging to T/Fb	Number of Leaf Cloud Points	Foliage Clump Volume (m3)/ Projection Area (m2)	Number of Leaves [29]	Leaf Area (m2)/LAI	Estimated Leaf Area Density (m2/m3)	Gap Fraction
				(Our Method/ Field Measurement)
Tree 1	A(Fb)	20832	29.28/3.26	1157/1274	8.74/2.68	0.30	0.42
	B(Fb)	17411	27.65/2.87	967/1027	7.30/2.54	0.26	0.53
	C(T)	21548	28.86/3.38	1197/1007	9.04/2.67	0.31	0.48
	D(Fb)	38216	49.12/4.23	2123/2242	16.03/3.78	0.33	0.43
	E(Fb)	19608	26.78/3.21	1089/1026	8.23/2.56	0.31	0.39
Tree 2	A(Fb)	11410	22.75/2.13	543/609	4.32/2.02	0.19	0.63
	B(Fb)	14821	20.34/2.30	706/789	5.62/2.44	0.28	0.61
	C(T)	8852	24.70/1.81	421/494	3.35/1.85	0.14	0.73
	D(Fb)	11884	25.05/2.23	565/487	4.49/2.01	0.18	0.75
	E(Fb)	41282	55.14/5.11	1966/2118	15.64/3.06	0.28	0.57

Note: T: trunk; Fb: first-order branch.

).

3.2. 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.

3.3. Analysing the Wind-Related Parameters in Forest Plots of Different Clones

According to the hurricane classification criteria in China, hurricanes are divided into six levels in accordance with the maximum wind velocity: Tropical depression (10.8~17.1 m/s); Tropical storm (17.2~24.4 m/s); Severe tropical storm (24.5~32.6 m/s); Typhoon (32.7~41.4 m/s); Severe typhoon (41.5~50.9 m/s); and Super Typhoon (>51.0 m/s). Our paper researched severe tropical storm Mirinae, which landed in Wanning, Hainan, on July 26, 2016. The study retrieved the aerodynamic parameters in the rubber tree forests of the Chinese Academy of Tropical Agricultural Sciences Proving Ground in Danzhou. When Mirinae landed, the maximum wind velocity near the hurricane eye was 25 m/s. It then moved northwestward and arrived in Danzhou at 6 a.m. the next day with the maximum wind velocity near the hurricane eye 18 m/s. Since the rubber tree experimental farm is approximately 15 km away from the hurricane eye, two forest plots composed of rubber trees of clone PR107 and CATAS7-20-59 were placed under the load of a hurricane with a speed of 15 m/s. We studied three trees in the second (middle) row and obtained the velocity, dynamic pressure and turbulent intensity results in the horizontal and vertical profiles, as shown in Figure 10 and Figure 11.

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.

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 (Figure 10c and Figure 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 prone to generate branch breakage and serious defoliation. However, in Forest Plot 2 (Figure 10f and Figure 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.

4. 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].

5. Conclusions

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

• 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.
• 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.
• 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.