Age Dynamics of Wind Risk and Tree Sway Characteristics in a Softwood Plantation

Quantifying wind risk in planted forests is critical to ensure the provision of wood-based biomaterials and forest ecosystem services in the future. Tree sway dynamics play a major role in determining their vulnerability to wind storms. While structural engineering theory can greatly assist in predicting tree and branch motion accurately, the effects of forest aging and silvicultural treatments on motion and the likelihood of mechanical failure are still poorly understood. In this study, we simulate the growth of an individual model tree and its architecture for a managed, even-aged pure Cryptomeria japonica D. Don. stand. The dynamic behavior of the tree is predicted for a rotation period at a 5-year interval using finite element analysis. Results show that tree natural sway frequency is subject to a 5-fold decrease over the tree's lifetime. The tree is thus more susceptible to absorb turbulent kinetic energy from the wind as it grows older. The damping ratio, which measures how rapidly the tree can return to its rest position after being deflected, is predicted to increase by a factor of four. By simulating the mechanical response to a wind gust, the model predicts a quick transition from a low risk of failure in the juvenile phase to a high-risk situation in the mature phase. The transition occurs between ages 15 and 30 when canopy closure takes place, when juvenile wood formation stops in the lower stem region and when the first thinning is carried out. Wind risk plateaus after the transition period. The age-risk relationship shown by the model is consistent with post-typhoon damage patterns documented by inventories in Japanese planted forests, albeit with a sharper transition.


INTRODUCTION
Trees are the longest-living organisms on Earth. To ensure that longevity, they must withstand the forces exerted on them by the physical environment throughout their lifetime. Although tree mechanical failure is a natural disturbance process important to forest dynamics and biodiversity (Thompson et al., 2009;Mitchell, 2013), large-scale damage is mostly detrimental in forestry because of the associated economic losses (Gardiner et al., 2013) and forest health issues (Stadelmann et al., 2014). As trees grow in size and exposure, they are subject to continuously increasing mechanical forces that should eventually compromise their likelihood of survival. Studies have shown that tree allometric patterns ensures that mechanical stability is less limiting for tree development than other constraints such as hydraulic transport (Koch et al., 2004;Niklas and Spatz, 2004). Often, those studies only consider elastic instability (buckling) under self-weight and neglect other agents and processes that can cause failure.
As top-loaded, tall and slender mechanical structures, forest trees are inherently unstable. To counteract this, they have evolved adaptive growth mechanisms limiting the effects of chronic and acute mechanical loading. A large focus in the study of tree biomechanics has been posture control and gravitropism, an active reorientation process to keep the tree upright thus minimizing the effects of weight (Fourcaud et al., 2003;Coutand et al., 2007;Alméras and Fournier, 2009). Other key adaptive mechanisms are thigmomorphogenesis (Jaffe, 1973) and thigmotropism (Telewski, 2012). They involve regulatory cambial growth to increase resistance and decrease exposure to mechanical loads. Trees also display passive mechanisms such as streamlining (Harder et al., 2004;Butler et al., 2012) and branch damping (Théckès et al., 2015). They allow reducing excess drag and improving kinetic energy dissipation in the swaying tree, respectively. They also make aerial architecture, (i.e., the geometrical and topological configuration of the branch system) a significant actor in tree biomechanical function and safety. The contribution of all those mechanisms is essential to ensure long-term tree stability.
Wind is a major abiotic agent of damage in forest plantations. While mechanical stability to wind loading has been studied in relation to growth (Ancelin et al., 2004;Gardiner et al., 2008), the effects of aging on tree vibrational properties and aerial architecture are rarely included. Instead, the dynamic behavior, which is essential to quantify the risk of windthrow and windbreak, is often investigated at a fixed point of a tree's lifetime (Moore and Maguire, 2008;Sellier and Fourcaud, 2009). As a result, the effects of canopy closure, silviculture, and juvenile wood formation on tree structural dynamics are poorly understood. That knowledge is critical to evaluate the vulnerability of forest plantations to extreme storm events. Extreme storm systems are expected to become regionally more intense and more frequent in future climate projections (Ulbrich et al., 2008;Gastineau and Soden, 2009). Because those events are rare, it is difficult to expand our knowledge about the damage they cause to forests.
In this study, we investigate how tree sway properties and the dynamic response to wind loading vary with tree age on a theoretical basis. The case study is an individual Cryptomeria japonica D. Don. tree representing managed, even-aged pure stand conditions. As Japan is annually subject to several tropical cyclonic storms (typhoons) and has ∼10 million hectares (30% of land area) of planted forests (MAFF, Japan, 2017), extensive information on tree growth and wind damage was available to support model development and evaluate the predicted age patterns of wind risk.

MATERIALS AND METHODS
The aerial system of a C. japonica tree was reconstructed from age 5 (sapling) to age 55. The common age of harvest in Japanese forestry is 60. The model tree was assembled using distinct experimental datasets and relationships for tree stem, architecture, foliage, and wood properties. A complete dataset was unavailable, even for a single age. The dynamic behavior of the tree was evaluated under impulse loading and wind loading using finite element analysis.

Tree Stem
The shape of the stem corresponds to a measured C. japonica tree named Sugi 2 (Minato et al., 1989). The tree was planted in the Wakayama experimental forest (135.40E 33.4N, Honshu island, Japan). It was 53 years old, 26.8 m tall and had a DBH of 28.5 cm when it was felled. The stem taper (diameter variation over height) was reconstructed from the internal ring structure at a 5-years interval. Stem taper at age 55 was obtained in this study by extrapolating ring width using the last known growth rates.

Crown Structure
The branching structure was modeled using the crown patterns measured in 9, 16, and 22-year-old stands (Hashimoto, 1991). Those ages represent the main phases of crown development from a young crown covering the full length of the stem to a mature crown after canopy closure. The number of branches borne per unit length of stem was given by a random Gaussian distribution (22.5 ± 3.3 branch m −1 ) (Hashimoto, 1991). Branches were uniformly positioned over the growth unit as C. japonica does not form whorls. Each branch was assigned a random growth potential at formation which acts as a scaling factor for initial length and growth rate. The upper crown length was defined as the topmost 4 m of the crown (Kajihara, 1976). In the lower crown region, branches were assigned a probability p = 0.1 to die every year because of limited light availability. Overall crown length increased with age up to a maximum of 9 m (Kajihara, 1976). The azimuthal orientation of branches was predicted using a phyllotactic pattern with a divergence angle equal to 139.5 • . The branch angle with the vertical is higher in older trees (Hashimoto, 1991), which was represented using θ = 71 − 0.5n where n is tree age at time of branch formation.
Branch primary growth was modeled using an age-dependent yearly length increment: L = a n b e c n + d (Hashimoto, 1991), with a = 0.06; b = 2; c = −0.5; d = 0.2. Branch basal diameter, D, was increased with branch length as D = 0.28L. The number of trees per hectare (stocking), s, ultimately limits the lateral extent of branches. For a circular crown, maximum length is equal to L max = 1 cos θ 10 4 π s . A thinning event causes s to decrease, allowing branch elongation to resume to take advantage of newly available space. At the Wakayama site, thinning events took place at age 20, 27, and 47 (Minato et al., 1989). In this approach, crown size expansion can resume or accelerate after thinning occurs. Modeling thinning was essential to synchronize crown growth with stem growth. The architecture of the model tree during growth is shown in Figure 1.

Foliage Properties
Tree's foliage mass, m F (n), is required to simulate tree mechanical behavior at any age n. It allows the inertia, the weight and the sail area of the crown to be calculated. In absence of repeated measurements during the growth of the tree Sugi 2, foliage mass was estimated using allometric models. Two models were considered. Both predicted foliage mass in relation to stem size or mass, which were measured (Minato et al., 1989). Only these two models provided a realistic trend of foliage biomass among considered models.
The first model relates directly foliage mass to stem diameter at breast height, D (cm). The model (Hashimoto and Gyokusen, 1995) predicts foliage mass as: It was retained as the best model describing the functional relationship between foliage mass and shoot cross-sectional area derived from plant hydraulic principles. The second model, purely empirical, was developed to obtain the most realistic allometric for the studied species. It predicts the foliage-tostem mass ratio as a function of age. It was obtained by fitting a power law on the biomass ratio using measurements across nine C. japonica forest sites at multiple times during the rotation. Those sites were Gunma/Ibaraki, Saitama, and Nara (Ando et al., 1968), Chiba (Negishi et al., 1988), Gifu (Watanabe and Moteki, 2007), Miyazaki, Nagasaki, Oita (Tadaki et al., 1965), and Toyama (Kato et al., 2005). They cover a broad range of silvicultural and environmental conditions, including different stocking densities and thinning schedules. The power law was fitted using Standardized Major Axis (SMA) regression as provided by the smatr R package (Warton et al., 2012). With both foliar models, the effects of thinning on total foliage mass result from the effects of thinning on stem growth (diameter or mass).
The amount of foliage borne by each branch was a fraction of the total crown foliage weighted by the basal diameter of the branch. Foliage was then distributed over the five most recent growth units of each branch to match the average lifespan of needles in C. japonica (5 years). The vertical distributions of foliage mass were in good agreement with profiles measured at ages 9, 16, and 22 (Hashimoto, 1991).
The foliage area, A F , is used for calculating the drag force applied to the crown of the tree. For an individual tree, it is approximated by A = m F (n) a F , where a F = 1.667 [m 2 kg −1 (green weight)] is the specific leaf area (SLA) (Tadaki, 1966).

Wood Properties
Three material properties were used to define wood mechanical behavior in the model. Those are the longitudinal Young's modulus, E (GPa), the density ρ (kg m −3 ), and the damping ratio, ξ 0 . E and ξ 0 are often known for dry wood. However, their values differ for green material in trees. Green wood is typically less stiff (lower E) and more viscous (higher ξ 0 ) than dry wood.
Most softwood species exhibit a typical radial pattern (TRP) of wood properties. The TRP consists of an increase, or a decrease depending on the property considered, followed by a plateau as the age of wood formation increases. In C. japonica, E follows the TRP (Nakatani et al., 1988;Zhu et al., 2005) and can be expressed as a function of the tree-ring number k (Zhu et al., 2005) as E k = c G (2.369 + 2.263 ln k ) where c G = 0.7 is the green-to-dry ratio E (arbitrary value for a 30% reduction in stiffness). Wood forming higher in the tree is stiffer than in the lower stem region (Yamashita et al., 2009). The effect of the height at which the wood is formed z was modeled by E k, z = E k (0.722 + 0.169 ln (5.204z) ). The cross-section of any stem or branch element is composed of several rings. The effective modulus of elasticity (MOE) of a cross-section was calculated as the sum of each ring's E value weighted by ring's moment of inertia: where R k is the outer radius of the kth ring. Unlike E, basic density remains relatively constant with cambial age in C. japonica (Minato et al., 1989;Yamashita et al., 2009). It was thus modeled as a homogeneous quantity. The green density value, 1,213 kg m −3 , was obtained using the average basic density (ρ= 438 kg m −3 ) and the average moisture content (w = 1.77) of Sugi 2 (Minato et al., 1989). The density value of woody elements bearing foliage was increased by foliage mass divided by the element's volume. The damping ratio of wood depends on the ratio of Young's modulus to basic density (Ono and Norimoto, 1983;Brémaud et al., 2013). For the kth ring, it was predicted using: where the specific density γ is equal to the ratio of basic density to water density. The factor d G = 1.27 corrects the material viscosity for the fact it is greater for saturated cell walls. The value of d G was calculated from damping experiments in relation to moisture content (Obataya et al., 1998). The effective damping ratio of a cross-section ξ 0 was calculated as the sum of each ring's ξ 0 value weighted by ring area. Material damping was then modeled using Rayleigh's model, which linearly relates damping to mass using the coefficient α R = 4πf 0 ξ 0 where f 0 is the natural sway frequency of the tree.

Numerical Model
The mechanical behavior of the tree was simulated using a finite element (FE) model developed and tested to study tree dynamics (Sellier et al., 2006(Sellier et al., , 2008. The model uses ABAQUS software with the AQUA module (SIMULIA, Dassault systèmes, France). Tree stem and primary branches are meshed using beam elements. Foliage was modeled by modifying the drag area and the density of elements bearing foliage. Dynamic analyses were carried out by direct integration of the equations of movement using an implicit solver. Unlike with a modal procedure, the geometry of the structure is updated at every time step of a simulation using non-linear analysis. This is essential to simulate large structural displacements, which may occur in response to gale force winds. Furthermore, geometrical non-linearity accounts for the streamlining phenomenon, (i.e., a reduced drag area when the crown is deformed by wind action). There is no need for ad hoc modifications of the drag equation.

Stability Criteria
Two types of simulations were carried out with the FE model. For both types, the mechanical behavior is simulated every 5 years. The first experiment was a virtual pull-and-release experiment. The tree was initially displaced by applying static wind loading on the tree's crown, then released suddenly to induce movement. A 100 tree architectures were randomly generated and tested. The crown features affected by random state were: number of branches borne by a growth unit, height of insertion, growth potential, and branch survival in the lower crown. The foliar model fitted on the experimental datasets was used to predict foliage biomass in this experiment because of its stronger empirical basis. Two important sway characteristics were calculated from the free oscillation response: • f 0 , the natural sway frequency of the tree. It is the inverse of the period required by a tree to sway back-and-forth once. f 0 is also referred to as the frequency of the fundamental (i.e., of lowest frequency) mode of vibration. Although the presence of modes of higher frequency has sometimes been observed in tree motion, the fundamental mode typically dominates the response. • ξ 0 , the damping ratio of the fundamental mode. ξ 0 measures the decrease of sway amplitude every cycle. For example, a structure would perpetually oscillate for ξ 0 = 0; the amplitude of lateral movement is divided by two every sway when ξ 0 ∼ 0.11; sway amplitude is divided by four every cycle when ξ 0 ∼ 0.22; and there is no oscillation (single cycle) when damping is critical (ξ 0 = 1). ξ 0 was calculated by fitting the theoretical motion of a simple damped oscillator to simulated tree response using least-squares optimization.
In the second numerical experiment, the behavior under wind loading was simulated. Because of the computational cost of those simulations, only one tree architecture was tested under wind loading: that of the tree with the sway frequency and the damping ratio closest to the mean values obtained in the first experiment. Simulations were done with both foliar models (section foliage properties). Two criteria were used to characterize the tree response.
• M 0 , the bending moment at the base of the tree. If M 0 exceeds the anchorage strength of the soil/root plate, structural failure by uprooting occurs. Anchorage strength is well-predicted by stem mass m S or, equivalently, volume (Moore and Quine, 2000;Cucchi et al., 2005;Kitagawa et al., 2010). The ratio M 0 m S was used to indicate the relative risk of uprooting. Other aspects not considered here, such as soil water logging or limiting rooting depth, may further reduce anchorage strength and thus increase the relative risk. • σ (z) the longitudinal stress along the tree stem. Breakage can occur at any point on the stem where σ exceeds the Modulus of Rupture (MOR). Thus, it is max (σ ) and not the mean stress that is relevant to analyze the risk of structural failure by windbreak. A single point of failure on the stem is sufficient to cause failure. MOR = 30 MPa is used as reference for green C. japonica wood. Natural variation is ∼± 50% of that value (Nakatani et al., 1988). It is consistent with MOR values from bending tests on green logs from a 33-year-old plantation (mean 34.9 MPa, min. 25.8 MPa max. 49 MPa, n = 36) (Suzuki and Ito, 1994).

Wind Loading
Wind exerts a force called aerodynamic drag (Cullen, 2005) on any tree element that is function of the speed of the wind relative Frontiers in Forests and Global Change | www.frontiersin.org to the structure: where ρ is the air density (1.25 kg m −3 ), c D is the drag coefficient (0.8), A is the projected area exposed to the moving fluid (m 2 ), u (z, t) the horizontal wind speed (m s −1 ) andq the speed of displacement (m s −1 ) of the tree element considered. A was calculated as the product of foliage area per unit length and element's length projected on the vertical plane normal to wind direction. The value of A was updated during simulations using the current element orientation. u (z, t) was decomposed into a spatial component u(z) and a temporal component u(t). To derive the spatial component of wind speed inside a forest stand, one must first define the wind speed when unperturbed by the forest (at the edge). Horizontal wind speed varies with height (Peltola et al., 1999;Gardiner et al., 2008). At the edge of a forest stand, the profile of u(z) is logarithmic (Ancelin et al., 2004): where h is the mean height of the canopy and z 0 is the roughness length (z 0 = 0.06h). u E h was calculated for each tree age using the above equation with u E (30) = 20 m s −1 for reference. At a distance more than 8-10 h away from the edge downwind, forest stand conditions apply (Dupont et al., 2011). Due to the presence of vegetation, the wind profile in the stand u S (z) is rapidly attenuated below the top of the canopy. Here the wind profile is represented by a power law: where u S h was obtained with a logarithmic profile above the canopy and u E = u S at 1 km altitude (Ancelin et al., 2004). The attenuation coefficient α was chosen as: where z C is the height at the base of the crown and the ratio u (z C )/u(h) is equal to 0.2. The resulting wind speed profiles are shown in Figure 1.
The time-dependent component of wind speed u (t) corresponds to a single gust event. In forest canopies, gusts are characterized by periodic ramp patterns of turbulent quantities with large deviation from the mean (Raupach et al., 1996). In the case of horizontal velocity, a gust is associated with air movement slowing down followed by a strong acceleration (Collineau and Brunet, 1993). In this study, corresponds to a 30-s event measured below canopy top (Collineau and Brunet, 1993). A maximum velocity equal to twice the mean is typical of wind storm conditions (Oliver and Mayhead, 1974).  Figure 2 shows the evolution of stem biomass for 60 years in nine different C. japonica plantations. At a young age (9-yearold), allocation was prioritized toward foliage production with nearly 60% of biomass allocated to it. This was measured at two forest sites: Saitama (Nishikawa management, 4,000-4,500 trees ha −1 ) and Gunma/Ibaraki (National Forest management, 3,000 trees ha −1 ) (Ando et al., 1968). The fitted model predicted up to 75% of biomass allocated to foliage for a 5-year-old tree. The situation changed rapidly: stem biomass was three times higher than foliage biomass at age of 20. The relative stem biomass kept increasing as the forest grew older, but more slowly. At harvest age, stem biomass was ∼6 times the foliage biomass.

Above-Ground Growth Allocation
The trend of allocation between foliage and stem was consistent for all sites despite varied initial stocking (up to 10,000 trees ha −1 at the Nara site under Yoshino management), cultivars and geographical locations (Kyushu and Honshu islands). Unlike the foliage:stem biomass ratio, the absolute foliage biomass productivity is highly variable across sites (not shown). The allometric relationship (Hashimoto and Gyokusen, 1995) predicted a evolution of allocation similar to the national trend obtained with the fitted model although it overestimated the amount of foliage at early age and underestimated it later on. It performed remarkably well for a simple relationship based on functional hydraulics without any age information. It also predicted more short-term variation in foliage biomass.

Natural Sway Frequency
The modeled tree displayed a general decrease in natural sway frequency with age ( Figure 3A). The decrease was observed whether the tree was debranched, defoliated, or complete. There was an exception for the complete tree which, at age 5, was found to sway more slowly (f 0 = 0.85 Hz) than at age 10 (f 0 = 0.93 Hz). Afterwards, f 0 decreased until it reached 0.19 Hz at age 55. Defoliation had for effect to increase f 0 with respect to the intact tree. The relative increase was greatest at age 5 (52%) and less but notable (19-30%) at older ages. Branch removal caused an additional increase in f 0 (i.e., the tree stem swayed faster without any crown). The increase due to debranching in addition to defoliating was smallest (12%) at age 10-15 and largest (27%) at age 35.

Damping Ratio
The damping ratio ξ 0 increased with tree age (Figure 3B). Damping was relatively low (ξ 0 = 0.05) when the tree was young (15-yo or less). It increased until the tree was 40-yo and stabilized at ξ 0 ∼ 0.2. The amount of variation in the damping ratio due to architectural differences increased as the tree got older. The quality of the fit used for determining ξ 0 was found unsatisfying at age 5 and at ages above 35 for different reasons. At age 5, a shift in natural sway frequency occurred while the tree was oscillating. A single degree-of-freedom (1-DOF) model cannot describe that phenomenon (Figure 4). At age 50, oscillation damping exhibited a distinctive non-linear behavior. Motion was strongly damped during the first cycle, but residual oscillations persisted for a long time afterwards. As a result, the 1-DOF model with linear damping under-then over-estimated the amplitude of tree sways during free vibrations. The r 2 value was high because the model performed well-quantitatively on average but less so qualitatively. This occurred from age 40 onwards. At age 20, a transient component was found in tree sway response (Figure 4). It probably resulted from the initial activation of a secondary vibration mode of higher frequency. The 1-DOF model could not represent that additional contribution and overestimated the amplitude of displacement during the first cycle. However, the quality of the fit is considered acceptable because the r 2 value was high and the signal perturbation was short-lived.

Sway Frequency and Crown Structure
Variability in branch positions in the crown caused variation in sway frequency. At age 5, f 0 ranged from 0.73 to 0.93 Hz FIGURE 4 | Relative lateral displacement of a fixed point on the tree stem during free oscillations at a young age (5), a mature age (20) and near harvest (50). Free oscillations were simulated by finite element analysis (blue line) for trees bearing individual first-order branches. A one-degree-of-freedom oscillator model (black line) was fitted to the simulated displacement time series to calculate the damping ratio. During the test, the initial deflection before release was applied to the tree crown using static wind loading.
depending on the crown arrangement, which corresponds to 22.6% of the mean. The range of variation was lower from age 10 onwards, between 2.5 and 8% of the mean. Figure 5 shows the relationship between f 0 (normalized by the sway frequency of the stem alone, f stem 0 ) and simple crown descriptors as the crown transitioned to a mature stage (age 25). At age 5, f 0 and crown mass m C were independent as f 0 varied by more than 20% for the same m C value (Figure 5A). At ages 10 and 15, no clear relationship was observed. At ages 20 and 25, a linear relationship was found. At those ages, it must be noted that f 0 was 30% less than for the stem alone whereas the crown represented 15% of tree mass. A better predictor of the influence of the crown on f 0 was found to be the change in height of center mass, z M , associated to the presence of the crown. The relationship between f 0 and z M were distinctively linear at all ages but 5 (Figure 5B). At 5, a crown of similar mass could either lower or heighten the center of mass depending on how the branches were positioned. As a result, f 0 values were quite dispersed. In older trees, a heavier crown always translated to a center of mass located higher than for the stem alone.

Response to Wind Loading
The bending moment and bending stress applied to the tree during a wind gust varied with tree age and size. The basal bending moment increased throughout the rotation period (Figure 6). The overall increase was three orders of magnitude for both the mean and the maximum values. Although the increase slowed down over time, the moment reached the highest value at the oldest age.
The magnitude of bending stress increased by a factor of 3-4 during the entire rotation (Figure 7). Unlike the bending moment, the mean bending stress peaked 5 years before the oldest age. The maximum bending stress caused by the gust event peaked even earlier, at age 40. Both mean and maximum stress values transitioned from low stress levels in the early phase to high stress levels that plateau later. In the case of the maximum stress, the transition was rapid, between age 10 and 25, whereas it took longer for the mean stress.
Both predicted basal bending moment and bending stress were higher when using the model describing the national trend of fractional foliage biomass. This stemed from the additional foliar biomass and crown area predicted by that model. Despite differences in load magnitude, the age dynamics of wind loading were similar for both foliar models. FIGURE 6 | Bending moment at the base of the tree induced by a wind gust as a function of tree age. The mean and maximum moment values recorded during the gust are given. Two model trees with different models of foliar allometry are considered: a fit of experimental allometric ratios and an empirical power law (Hashimoto and Gyokusen, 1995). Figure 8 shows peak tree deflection during a wind gust at three different tree ages. When the tree was 15-or 25-yo, the lateral displacement was moderate. However, the stress profiles were completely different. Bending stress was uniformly distributed along the length of the stem at age 15 whereas it is concentrated in the upper crown region at age 25. At age 45, the windinduced lateral displacement became large and the curvature more pronounced. Bending stress was high over most of stem length excluding the lower region extending from tree base and the very top of the tree. Figure 9 shows how wind risk evolved with the age of the plantation forest using the maximum loads obtained under wind loading (Figures 6, 7). Risk of structural failure is given as a ratio of load to notional resistances to windthrow and windbreak. The risk of windthrow ( Figure 9A) increased up to age 30 with both foliar models considered. The likelihood of overturning remained high in older trees without increasing further. The fitted foliar model yielded a risk of windthrow that was 10.9% higher on average, all other factors being equal, which is consistent with higher foliage area.

Age-Risk Relationship
The stress maximum in the stem did not reach the lower end of the MOR range in young trees ( Figure 9B). Because compression strength is between 50 and 60% of the MOR value (Chiba, 2000;Zhu et al., 2005), local compression failures could happen. After the age of 20, the maximum bending stress reached a level which implies windbreak would occur in trees with weak wood. In the later growth phases (age ≥ 30), the maximum stress reached a FIGURE 7 | Bending stress induced by a wind gust in the tree stem as a function of tree age. The mean and maximum stress values recorded during the gust are given. The mean value is averaged in time and along stem length. The maximum value is the overall maximum recorded at any stem height at any time during the gust. Two model trees with different models of foliar allometry are considered: a fit of experimental allometric ratios and an empirical power law (Hashimoto and Gyokusen, 1995).
value sufficient to induce stem failure in most trees. However, tree stems with the strongest wood and no point of weakness could theoretically resist those stress levels. Like the risk of windthrow, the risk of windbreak was found high and stable after a 10-year period of transition (age 15-25) from a low risk level. The same evolution of windbreak risk was obtained independently of the foliar model, except for an average increase in maximum stress of 10.6% associated to using the fitted model for foliage allocation.

Evolution of Natural Sway Frequency With Age
The model predicted a 6-fold decrease of natural sway frequency over the lifetime of the studied C. japonica tree. The overall range 0.18-0.94 Hz, lies within the global 0.16-1.39 Hz range known for trees (Moore and Maguire, 2004). Albeit shifted toward higher values, the predicted natural frequency evolved with age similarly to measurements: 0.62-0.78 Hz in 4-yo Pinus pinaster Ait. (Sellier and Fourcaud, 2005), 0.65 Hz in a 6-yo C. japonica (Suzuki, 2012), 0.18-19 Hz for 36-yo P. pinaster (Sellier et al., 2008) and 0.11 Hz in a 63-yo C. japonica (Suzuki, 2020). Age is useful to define a main trend in sway frequency evolution during growth. However, frequency depends principally on size. Silviculture through thinning schedules and stocking densities FIGURE 8 | Tree maximum deflection during a gust event at ages 15, 25, and 45. Grayed silhouettes represent trees at rest. Color scale, from blue (low) to red (high) corresponds to the relative magnitude of longitudinal stress under wind loading.
will modify the age-size relationship and, therefore, the agefrequency relationship.
The decrease of f 0 with age is consistent with mechanical theory which predicts that f 0 varies as the square root of the stiffness to mass ratio [i.e., a specific stiffness (∼ E 1/2 D L 2 for a cantilever beam)]. As the structure grows heavier and taller, mass scales faster than bending stiffness and the frequency decreases. Those effects are mostly size effects dictated by growth rates. Trees of different ages with similar dimensions will oscillate at similar natural frequencies. Age might have a direct effect on f 0 through the modulus of elasticity, E. The progressive stiffening of wood with age of formation would increase f 0 . That effect must be minor because it is insufficient to counteract the general decrease of f 0 observed in the numerical simulations.
The sway frequency helps in estimating the stability of an individual tree. Movement is amplified when the frequency of a dynamic load is similar the frequency of one mode of vibration of the structure (Wood, 1995). Trees are more sensitive to the turbulent kinetic energy of wind in the frequency range close to their natural frequency. In non-storm conditions, the peak frequency of forest canopy turbulence is less (0.02-0.05 Hz) than that of tree sways (Amiro, 1990). Higher peak frequency values equal to 0.1 Hz have been recorded during a typhoon (Yamanoi et al., 2005). Trees with the lowest f 0 values, (i.e., the older and larger ones), are more at risk of experiencing large sway displacements due to dynamic amplification and resonance. The decrease of f 0 with age is asymptotic (0.25 Hz at age 35 vs. 0.18 Hz at age 55). The transition to a less stable state, with respect to sway frequency, occurs rather early and rapidly. The windthrow risk indicator is given as the ratio of the basal bending moment to stem mass, which predicts well tree overturning strength. (B) The relative windbreak risk is given as the ratio of maximum mechanical stress in tree stem to the wood's Modulus of Rupture (MoR). The grayed area in (B) indicates ±50% of reference MoR; it covers most of the natural range of variation for green C. japonica (Nakatani et al., 1988). Results shown for both foliar models.

Evolution of Damping Ratio With Age
Between ages 5 and 40, the tree structure transitioned from being lightly damped (ξ 0 ≈ 0.05) to being heavily dampened (ξ 0 ≥ 0.2). That transition corresponds to a positive effect of aging on tree stability to wind forces. A higher damping limits the dynamic amplification of tree sways under dynamic loading more effectively. This positive effect contrasts with the negative effect of aging on frequency and wind load magnitudes. Values predicted for ξ 0 lie within the 0.03-0.25 range documented by (Moore and Maguire, 2004). However, measured ξ 0 values are more commonly found between 0.05 and 0.1. Furthermore, no relationship has been documented between ξ 0 and tree size and, a fortiori, age (Moore and Maguire, 2004). ξ 0 can also vary substantially among trees of similar age and dimensions [e.g., 0.05-0.12 in pine saplings (Sellier and Fourcaud, 2005)]. Because of lack of data, it is not possible to assert that the simulated increase of ξ 0 occurs in nature.
The predicted increase of damping with age is consistent under the hypothesis that aerodynamic drag is the main source of damping in trees. The exponential increase of foliage biomass in the early stage of tree growth would translate to a rapidly increasing aerodynamic damping. That increase would cease after canopy closure unless crown growth is released by thinning events. Once crown size and foliage mass have stabilized after release, the aerodynamic damping would plateau as it is observed in the numerical simulations. The higher ξ 0 value predicted at age 50 is likely a response to thinning at age 47. Trees also dampen motion thanks to branch vibrations (Sellier et al., 2006;Spatz et al., 2007). The efficiency of branch damping does not depend on crown size but on a fine-tuning of branch natural frequencies (Scannell, 1984;Théckès et al., 2015). Viscous damping, which dissipates motion though internal friction in wood, is known to have a minor contribution to overall damping (Milne, 1991;Sellier and Fourcaud, 2005). Viscosity should decrease for wood formed later in tree lifetime declines because of the inverse relationship between viscous damping and specific modulus (Ono and Norimoto, 1983).
Two important aspects have emerged about damping in this study. Firstly, the absolute damping rate varied little with tree age. It means that the tree deflected by a wind gust remains swaying for a similar time window (20-30 s) whether it was 5, 20, or 50 years old. During that window before the tree returns to rest, it is temporarily more exposed to the arrival of another gust. Secondly, the nature of tree sways changed as the tree grew. At a young age, tree's response exhibited a frequency shift. At an intermediate age, a secondary mode of vibration was excited. That mode limited the lateral displacement during the first cycles. At an old age, damping appeared quadratic (i.e., drag-dominated) with a relatively higher dissipation of initial cycles followed by persistent small-scale vibrations. In all those cases, the notional simple oscillator model used for calculating ξ 0 was insufficient to represent the aspects of damping at each age. However, ξ 0 remains useful to provide a gross estimate of the structure's ability to dissipate motion.

The Role of Crown Architecture in Tree Dynamics
From the numerical experiments manipulating the crown status, we find that the presence of branches and foliage slows down tree sways. In absence of any crown, at age 5, the tree swayed twice as fast as the complete tree. A similar ratio was measured for a 4-years-old P. pinaster sapling (Sellier and Fourcaud, 2005). The frequency ratio between the debranched tree and the complete tree continuously decreased during growth until it was 1.25 at harvest age. Thus, the crown's influence on tree sway frequency is reduced with tree age. The crown slows down tree motion because it adds mass onto the stem without contributing any stiffness. Despite that, the relationship between crown mass and frequency is weak. Most of the crown must be removed to significantly increase f 0 (Moore and Maguire, 2005). It is also found that the relationship between crown and sway frequency was better represented when the spatial distribution of crown mass was taken into account, for instance by using the height of the center of mass as a predictor of f 0 changes. The effects of branch distribution along the stem on natural sway frequency was found to be age-related too, with the youngest trees affected most. Integrated silviculture could become a tool for wind damage risk management by influencing the sway frequency of planted trees because it can change crown structure and allocation between compartments (Ando et al., 1968).
Even the center of mass was a poor predictor of the influence of crown on f 0 at age 5. At that age, the spatial distribution and the amount of crown biomass were loosely coupled. All other things being equal, a sapling with a heavy crown could have a center of mass at the same height as a sapling with a light crown. Inversely, two saplings with crowns of equal mass could have centers of mass at a different height. It is possible to accurately predict f 0 in the case of saplings, but those predictions require to describe tree's aerial architecture extensively (Sellier et al., 2006). This is in extenso also the case of mature opengrown angiosperm trees, where large branches have a significant impact on the dynamics of the whole tree (Rudnicki et al., 2012;Ciftci et al., 2013). The emergence of laser scanning technology to measure tree architecture will play a major role to systematically perform complex analysis like finite element modeling (Jackson et al., 2019) and to characterize the dynamical properties of trees during their development.

Aerodynamic Behavior: Juvenile and Mature
Wind forces continuously increase as the tree grows larger and taller. The height growth translates to exposure to higher wind speed and an increased lever arm for the bending moment. Foliage growth and branch growth cause greater sail area and drag. Despite tree dimensions and exposure continually increasing because of growth, the increase in wind forces becomes gradually smaller in the late stages of the rotation. This behavior may result from the shift in carbon allocation from foliage to woody biomass due to aging. The shifted allocation translates to relatively less crown sail area and a stronger stem. If that is the case, then there is a clear need to capture precisely allocation trends and to investigate the role of forest management on those trends in future studies of wind risk. We have presented two allometric models that predict reasonably well foliage mass and crown area across a wide range of sites, even for a single individual tree. They are a practical alternative to repetitive foliage measurements for a plantation lifetime. Hashimoto and Gyokusen (1995)'s model is interesting in that it describes the functional relationship between foliage area and sapwood area promoted by the pipe model theory (Lehnebach et al., 2018). In the future, conceptual models of allometry based on the constant-stress theory (Morgan and Cannell, 1994) or source-sink relationships (White et al., 2016) could be used or tested.
In the juvenile tree, mechanical stresses are approximately uniform along the length of the stem. This agrees with the uniform stress theory that is a well-known paradigm in plant biomechanics. Stresses induced by wind loading are expected to be uniform below the crown and above the butt swell (Morgan and Cannell, 1994). The mature tree also presents a long stem section where bending stresses are uniform. This corresponds to the stem region where breakage is expected to occur (Chiba, 2000). That section extends well into the lower crown region. The modulus of rupture of wood increases with height in trees (Nakatani et al., 1988), which might be a beneficial trait to withstand higher stress levels in the stem region where branches are borne. Breakage of the top part of old trees during typhoons sometime occur in old C. japonica specimens growing in natural forests (Takashima et al., 2009). That type of windbreak could be associated to a decurrent tree form, however, and is atypical in plantations. The tree at age 25 presented a significant deviation from uniformity with very high bending stresses near tree top. The exception can be explained as a response to thinning at age 20 with crown growth in width temporarily exceeding stem radial growth.

Age-Risk Pattern and Consequences for Future Forestry
The risk of wind damage was found to strongly depend on which age class the tree belongs to. In the simulations, trees <20-years-old are submitted to low-magnitude wind loads. All the criteria considered, except for damping, indicate that juvenile trees can be considered as stable under wind loading. During an intermediate phase (age 20-30), the tree transitioned from a situation of low risk to one of high risk. Several events occurred during that period: LAI peaked, the canopy closed, wood properties plateaued, and two thinning events were applied. Any of those events would be detrimental to stability on their own. Because all those changes were concurrent, it could explain why the transition to a situation of high risk was particularly rapid and pronounced. In future studies, it would be interesting to quantify how successive thinning events interact with canopy growth and, by releasing growth, if they contribute to an artificially heightened or precocious risk of wind damage. This would allow the onset and magnitude of high wind risk to be tied with silvicultural planning and productivity curves. This and the role of planting density could be investigated using a combination of growth and yield models and theoretical allometric models.
In the second half of the rotation period (after age 30), the predicted level of wind risk remained approximately constant despite the continuously increasing wind loading on the structure. It appears that the strength of the stem and the strength of the root system increased in proportion to the applied loads. Dedicated models of stem and root plate mechanics should be employed to investigate the scaling of strength with age in more detail.
Wind damage recorded in planted forests after typhoon events has a distinctive pattern. Minimal damage is reported for trees 20years-old or younger. Later, the amount and likelihood of damage increase to a peak value followed by a decrease. Numerous damage inventories reported this "bell-shaped" pattern (Nakao et al., 1993;Fujimori, 1995;Abe, 2005;Sato and Abe, 2006). It is observed across different regions, from Kyushu to Hokkaido, and across different planted species [C. japonica, C. obtusa, Larix kaempferi (Lamb.) Carr.]. The age of peak damage can vary significantly across sites from 35 (Nakao et al., 1993) to 80 (Sato and Abe, 2006), commonly between 50 and 60 (Fujimori, 1995;Abe, 2005). Deviation from the general pattern may happen. For example, important damage has been recorded in the sub-20 age classes at one site (Nakao et al., 1993). Damage inventories for natural forests have shown an altogether different age-risk relationship with a steadily increasing likelihood of damage over a 100 of year (Fujimori, 1995). Similarly, in the natural forests of the Pacific northwest (North America), wind becomes a significant cause of tree mortality only after age 20 and becomes the primary cause in age classes older than 100 (Franklin et al., 1987). The differences between natural forests and plantation forests indicate that the age-risk relationship shown in this study is not inherent to tree's increase in size during growth in general. Whether the influence of forest management on risk derives from the impact on tree growth (slenderness, taper, allocation) or the impact on forest structure heterogeneity remains to be investigated.
In Japan, the dominant age class in plantation forests was 50-55 in 2017. It is of paramount importance to determine whether the age of peak damage likelihood has been reached to decide if the rotation period can be extended without increasing the vulnerability of forests to wind storms nationwide. Age 50 is approximately the mean reported age for peak damage but, because that age is variable, it is difficult to estimate if the risk associated to windstorms has already peaked. Although the present model can, to a large extent, reproduce the documented age-risk pattern, the lack of available growth data unfortunately did not allow predicting tree mechanical behavior beyond age 55. Assessing wind risk beyond the current traditional age of harvest remains difficult because experimental data on biomass spatial patterns (including architecture) is scarce past that age and because of the strong coupling between those patterns and tree aerodynamic behavior. Linking tree mechanical behavior to the age dynamics of forest aboveground biomass allocation as a function of climate and forest management is an important challenge for future research on wind-tree interactions.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
DS developed and implemented the growth model, performed the finite element analyses, analyzed the results, and wrote the manuscript. SS sourced published and unpublished experimental data sets, translated scientific literature, analyzed the results, and edited the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
The Japanese Society for the Promotion of Science (Grant no. S12707) provided a short-term fellowship funding to support DS 5-weeks stay at FFPRI, Tsukuba, Japan in 2012.