We developed two models to predict the probability of wind damage based on a hierarchical Bayesian method. The first included wind speed, wind direction, season, and mean stand volume as potential factors influencing wind-damage probability (model 1), and the second included wind speed, wind direction, season, and decade (model 2). Due to imperfect record collection, the values for decade-specific mean stand volume across entire UTCBF in model 1 were estimated, as were the effects of mean stand volume and other factors.

Model 2 did not address explicit factors representing the developmental or structural status of forest stands, but included a temporally autocorrelated factor (decade) that could indicate temporal trends in wind-damage probability. Details of the two models are given below.

The main assumptions (the probability distribution of damage events and the relationships between wind-damage probability and influencing factors) for model 1 are as follows:

D _i ~ Bern ( p _i ) (1)

logit( p _i ) = b ₀ + b ₁ ve _i + b ₂ vo ( de _i ) + e ₁ ( di _i ) + e ₂ ( s _i ), (2)

where D _i is a binary variable taking a value of 1 if wind damage occurred in the target stands during the i th wind event and 0 otherwise, and is assumed to follow a Bernoulli distribution Bern( p _i ); p _i is the probability that the i th wind event will cause damage in the target stands; conditions during a wind event are described by wind speed ( ve _i ; continuous variable expressed in m s ^-1 ), mean stand volume ( vo _i ; continuous variable expressed in m ³ ha ^-1 ) for a decade ( de _i ; discrete variable), wind direction ( di _i ; discrete variable), and season ( s _i ; discrete variable); e ₁ - e ₂ are functions describing the effects of wind direction and season, respectively, whose values are autocorrelated with their neighboring values (explained below in detail); b ₀ is an intercept term; b ₁ and b ₂ are the model parameters relating wind speed and mean stand volume to wind-damage probability, respectively; and logit is a logit function. Wind direction was included in the model because wind direction can modify the impact of winds depending on topographic and geographic conditions or the degree of acclimation of trees to wind ( Ennos 1997 ; Nilsson et al. , 2004 ). Season was included because it modifies the impact of wind through changes in soil conditions and pheno­logy, i.e. the presence or absence of deciduous leaves ( Nilsson et al. , 2004 ; Usbeck et al. , 2010 ; Brázdil et al. , 2018 ). Mean stand volume was included because the impacts of winds can vary with developmental status ( Peterson, 2000 ; Albrecht et al. , 2012 ).

We assumed that the values of e ₁ and/or e ₂ would be autocorrelated with those of neighboring levels ( e.g. e ₁ (NE) would be correlated with e ₁ (NNE) and e ₁ (ENE)), and we therefore used the intrinsic Gaussian conditional autoregressive model (intrinsic CAR model; Besag et al. , 1991 ) to describe these effects. The arrangements of levels for these variables were one-dimensional and circular; each level had two neighboring levels. Thus, the conditional probability distribution of the autocorrelated effect is expressed as:

To provide an example of applying the models described above, we reanalyzed longitudinal wind-damage data ( Abrams et al. , 2017 ) recorded at UTCBF (Fig. 1). UTCBF, which was established in 1894, was the first university forest established in Japan ( Research Section of Tokyo University Forests and Tokyo University Forest in Chiba, 1974 ). UTCBF comprises 2,226 ha of various forest types in the warm-temperate zone. Forest types included plantations of Cryptomeria japonica (L.f.) D.Don and Chamaecyparis obtusa (Siebold et Zucc.) Endl. and natural forests dominated by Abies firma Siebold et Zucc., Tsuga sieboldii Carrière, and evergreen broad-leaved trees. The topography is characterized as steeply dissected, with a dominant ridge/valley direction of north-south and steep slopes (>30 degrees) with primarily east and west aspects (Fig. S1 [suppl.]).

Tree damage caused by extreme weather events, primarily strong wind events such as typhoons but also snow and ice storms, were recorded and compiled by Negisi (1997) . That report covered a long time period (1900 to 1985), but lacked detailed information such as damage locations, consistent descriptions of damaged forests ( e.g. stand age, structure, or species), and quantitative damage measures. Only imperfect historical data describing developmental status ( i.e. mean stand volume) are available for UTCBF.

We used existing wind data for extreme weather events and disasters ( Choshi Local Meteorological Observatory, 1969 , 1987 ) recorded at Katsuura, the nearest weather station to UTCBF (Fig. 1). The daily maximum wind speed and direction for days with relatively strong wind (maximum wind speed >10 m s ^-1 ) were used as explanatory variables in the model. Although meteorological data recorded at Katsuura from 1961 onward are available online ( http://www.data.jma.go.jp/gmd/risk/obsdl/index.php ), wind data collected prior to 1961 are only available as hard copy records. For the post-1961 dataset, the frequency distribution of the strongest winds within each month had a distinct peak at SSW; 90% of the monthly maximum winds blew from S-WSW (Fig. S2 [suppl.]). Additional, long-term meteorological data recorded at four meteorological stations (Choshi, Osaka, Miyakejima, and Kobe, which are long-term stations; Fig. S3a [suppl.]) on Honshu (the biggest island in Japanese archipelago) and a small island near Honshu, and those recorded at meteorological stations near UTCBF were summarized to understand temporal changes in wind conditions (specifically, the number of days with a maximum wind seed >10 m s ^-1 ). At all long-term stations excluding Osaka, where wind speeds are generally low, presumably due to its inland location, the number of days with strong wind increased from 1900 to the period of 1945-1960 and then decreased from this period to 1980, with some fluctuations (Fig. S3b [suppl.]). For time periods later than 1961, data from meteorological stations in the proximity of UTCBF were included (Katsuura and Chiba; Fig. S3b [suppl.]). The observed change in wind conditions among distant meteorological stations, at maximum >500 km apart, suggests that this pattern was unlikely the result of local conditions around stations ( e.g. forest development).

We categorized wind direction ( dii ) into 16 levels (N, NNE, NE, ENE,… NNW), season ( si ) into 12 levels (closely corresponding to months), and decade ( dei ) into 9 levels (1900-1909, 1910-1919, …, and 1980-1989). We assumed non-informative prior distributions for almost all model parameters. Sampling from the posterior distribution of all parameters was performed using Markov chain Monte Carlo (MCMC) simulations with Stan software ( Carpenter et al. , 2017 ) called from R ( R Core Team, 2019 ) using the library RStan ( Stan Development Team, 2018 ). We followed Joseph (2016) in implementing the intrinsic CAR model in our Stan code. Stan codes for the wind-damage models are provided in the supplemental material (Appendices A1 and A2). Four independent MCMC chains were run, and 200,000 samples were recorded after a burn-in of 2,000 for model 1; 40,000 samples sufficed for model 2. The chains were thinned every 50 runs, yielding independent samples from a posterior of 16,000 for model 1. Ȓ values were used to assess convergence; Ȓ values <1.1 indicated that the four chains had converged. For each model parameter, a mean value was calculated from the sampled values and used as the point estimate of the parameter.

To assess whether the resultant models reproduced the major temporal pattern of wind damage, we cal­culated the wind-damage probability using records of extreme weather obtained from Katsuura. We summed the probability from four intervals with a fixed length (20 years) to obtain the expected number of days with wind-damage events in UTCBF for those intervals.

To demonstrate the utility of our wind-damage probability model in understanding forest history, tree regeneration information based on tree-ring data obtained from an old-growth fir-hemlock forest in UTCBF was used to compare predicted temporal patterns of wind-damage probability and observed tree regeneration. We compared the observed tree regeneration with the predicted wind-damage probability, rather than directly with the wind-damage record, because the damage records ( Negisi, 1997 ) were likely imperfect, as they may have been biased toward plantations and included wartime periods. A brief history of the target fir-hemlock forest is as follows ( Abrams et al. , 2017 ): Prior to the 1920s, the target forest had a long history of management using the coppice-with-standards me­thod ( Honda, 1912 , 1927 ; Research Section of Tokyo University Forests and Tokyo University Forest in Chiba, 1974 ). Under this method, shade-tolerant evergreen broad-leaved trees under the canopy layer were frequently cut (once every 20-30 years; Honda, 1927 ), creating opportunities for less shade-tolerant conifers to regenerate. Some regenerated conifers (standards) were reserved for timber production and to form a sparse canopy layer. Thus, several old conifers have been retained to the present day. The last clear-cut of broad-leaved trees was conducted in 1919. A detailed description of this forest is given in Abrams et al. (2017) .

In July and September 2013, twenty point locations (35°10'51-55"N, 140°9'11-38"E), located at appro­ximately 35-m intervals on a ridge in the forest interior, were used for dendrochronological sampling (Fig. 1). Across all twenty locations, we obtained 64 cores from five major species ( Abies firma, Castanopsis sieboldii (Makino) Hatus. ex T.Yamaz. et Mashiba subsp. sieboldii , Cinnamomum yabunikkei H.Ohba, Quercus acuta Thunb., and Tsuga sieboldii ) over a range of diameter classes. Two of the dominant conifer species ( Abies firma and Tsuga sieboldii ) are less shade tolerant, whereas three sub-dominant evergreen broad-leaved species ( Castanopsis sieboldii, Cinnamomum yabunikkei , and Quercus acuta ) are shade tolerant ( Abrams et al. , 2017 ). Trees were cored (one core per tree) at 1.37 m from the ground for age determination. The regeneration year of an individual tree was defined as the year wherein the individual reached 1.37 m in height. We compared the temporal distribution of re­generation in the major tree species at the study site and the expected number of wind-damage events to assess whether regeneration had occurred following major wind disturbance events ( e.g. typhoons).

All Ȓ values for the parameters included in models 1 and 2 were <1.1, indicating that all four MCMC chains converged in each model fitting. The posterior means and 95% Bayesian credible intervals for model parameters are provided Tables S1 and S2 [suppl.]. The results of model 1 are highlighted below because the widely applicable information criterion (WAIC; Watanabe, 2010 ) was nearly identical (105.2) for both models, and past mean stand volumes were estimated in model 1.

The predicted wind-damage probability for UTCBF was mainly determined by wind speed, wind direction, and mean stand volume. The probability increased with increasing wind speed (Fig. 2a), reaching a maximum at N-ENE directions and a minimum at SW (Fig. 2b). Wind-damage probability generally increased with increasing mean stand volume (Fig. 2c; 95% Bayesian credible intervals for b ₂ included 0). Season did not have a perceptible effect on wind-damage probability (Fig. 2d). The estimated mean stand volume generally increased throughout the study period, 1901-1990 (Fig. 3). The increasing pattern observed in the estimated mean stand volume was similar to that of the temporally autocorrelated effect of decade in model 2 ( e ₃ ( de _i ); Fig. S4c [suppl.]). The Bayesian credible interval for wind-damage probability was wider for model 1 than model 2 (Fig. 2 and Fig. S4 [suppl.]).

Temporal trends in the expected number of days with wind damage corresponded well with the number of days with recorded wind damage at UTCBF (Fig. 4). Strong winds that would likely cause wind damage were most frequent (12.4 days in a 20-year period) from 1946 to 1965, and least frequent (3.1-3.3 days in a 20-year period) before 1945. Extremely strong winds that were highly likely to cause wind damage (0.75 < p < 1.00) occurred between 1926 and 1965, but did not occur in the preceding period. Very similar results were obtained with model 2 (Fig. S5 [suppl.]).

Remarkable differences were found in the temporal patterns of regeneration between the two dominant coniferous species ( Abies firma and Tsuga sieboldii ) and the two subdominant broad-leaved species ( Cin­namomum yabunikkei and Quercus acuta ; Fig. 5). Tree regeneration peaked during 1886-1905 for the two dominant coniferous species. The regeneration of Abies firma had a second, lower peak during 1946-1965. The long-term temporal regeneration pattern and the expected number of wind-damage events differed considerably for the dominant conifers. Tree regeneration peaked during 1946-1965 for the two subdominant broad-leaved species ( Cinnamomum yabunikkei and Quercus acuta ). The regeneration for these species peaked during the periods where the expected number of wind-damage events was highest (Figs. 4 and 5). Another broad-leaved species ( Castanopsis sieboldii ) had two regeneration peaks, during 1926-1945 and 1966-1985.

The statistical model presented in this study reproduced the temporal patterns of past wind damage and mean stand volume in UTCBF (Fig. 3), and cla­rified relevant factors determining the probability of wind-damage.

We introduced temporally or directionally auto­correlated terms ( e ₁ , e ₂ , and vo in eqn. 2 for model 1; e ₁ - e ₃ in eqn. 5 for model 2). For these terms, possible problems caused by autocorrelation ( e.g. inflated sample size and underestimated variance), which exist in most field data, were avoided using our statistical approach. Even where longitudinal data of potential influencing factors ( e.g. mean stand volume in model 1) are incomplete, the missing information can be estimated by assuming that the factors are temporally autocorrelated (Fig. 2c). Moreover, temporally autocorrelated random effects ( e.g. e ₃ ( de _i ) in model 2) may account for a temporal pattern associated with change in unobserved factors (Fig. S4c [suppl.]). However, caution must be exercised in estimating unobserved variables when data are lacking. Appropriate assumptions are necessary for valid estimation, and estimating a large number of values often leads to wide credible intervals. Our proposed model framework has another advantage in its expandability; if additional data for factors influencing wind-damage probability are available, these models can be easily expanded by adding terms to represent these factors.

One of the shortcomings of statistical approaches to wind-damage modeling is the low applicability of the resultant model to situations that differ from those where the model was developed. Wind-damage probability models developed with longitudinal data in one area may not be useful for other areas. However, these models are expected to accurately predict the probability of future wind damage within the study area if the longitudinal data contain relevant patterns of wind events and these patterns do not change in the future. Therefore, models such as those presented here can be used in local forest management by providing future probabilities of wind damage. Further, they can provide vital information for assessing forest histories.

As expected, the probability of wind damage in UTCBF increased with increasing wind speed (Fig. 2). Similar patterns have been observed in many previous studies ( Peterson, 2000 ; Ulanova 2000 ; Usbeck et al. , 2000; Donis et al. , 2018 ). This is logical, as the stress (force/unit area) imposed on the canopy by wind is considered a function of wind speed ( Gardiner et al. , 2000 ).

Interestingly, we found that the probability of damage was higher with N-E than with S-W winds. This pattern may be attributable to the combination of topography, wind direction, and acclimation of trees to wind ( Ennos, 1997 ; Nilsson et al. , 2004 ). At UTCBF, ridges and valleys tend to run north-south (Fig. S1 [suppl.]). This may increase wind speed in valleys when winds blow parallel to the valley orientation ( Ruel et al. , 1998 ). Furthermore, the majority of strong winds blow from the south in this region (Fig. S2 [suppl.]). Therefore, trees are likely to have acclimated to southern winds, and may be vulnerable to wind from other directions, especially from the north, as wind speed may increase when winds blow parallel to the dominant valley direction ( i.e. from the north).

Season did not have perceptible effect on wind-damage probability, although seasons can alter the wind-damage probability for deciduous species, as the force received by the tree from wind is dependent on their leaf phenology (green or fallen; Nilsson et al. , 2004 ). However, forests in UTCBF are dominated by evergreen species, which likely explains why season was not an important factor for this study area.

Although the records of mean stand volume were incomplete and potentially inaccurate, model 1 es­timated the temporal change in mean stand volume and detected a positive effect of mean stand volume on wind-damage probability (Fig. 2c and Fig. 3). Broadly, the forests at UTCBF may have become more vulnerable to wind damage as trees grew taller and were therefore more exposed ( Albrecht et al. , 2012 ; Donis et al. , 2018 ). Although no information on forest development stage was used in model 2, the temporally autocorrelated random effect ( e ₃ ) detected a temporal pattern similar to that of mean stand volume (Figs. 3 and S4c). This indicates that the average developmental status across UTCBF may be the major driving factor in the observed temporal pattern of wind-damage probability, which was not explained by wind speed or direction. We therefore suggest that the terms representing dynamic forest developmental status be included in future models that aim to analyze longitudinal data.

If tree regeneration and mortality rates are cons­tant in a population, the frequency distribution of regeneration year (a mirror image of age structure) should become an exponential function. However, temporal variability in regeneration or mortality can change this distribution from the exponential function. For example, if regeneration peaks in a given period, age structure shows a corresponding peak. Among the frequency distributions of tree regeneration year for major species in the target fir-hemlock forest, none was similar to an exponential function (Fig. 4), indicating that regeneration and/or mortality were not constant.

A remarkable temporal pattern was observed in the frequency distributions of tree regeneration year for the two dominant conifers ( Abies firma and Tsuga sieboldii ), which peaked during 1896-1905 (Fig. 4). This pattern differed not only from the exponential function but also from the temporal pattern of the expected number of wind-damage events (Fig. 3). As described in Abrams et al. (2017) , this regeneration pattern was likely created by historical human management and not by wind disturbances. The majority of fir-hemlock forests at UTCBF had been managed using a coppice-with-standards method until the 1920s. In this method, conifers in the canopy layer (standards) were reserved for timber production, mostly for large-scale construction, while evergreen broad-leaved trees under the canopy layer were used for fuel wood and charcoal as short-rotation sprout coppice (cut once every 20-30 years; Honda, 1927 ). This relatively frequent cutting of understory broad-leaved trees gave less shade-tolerant conifers the opportunity to regenerate, and once regenerated, conifers can be long lived as reserved standards. The regeneration year of existing conifers in the target fir-hemlock forest was dated as far back as 1810 and peaked during 1886-1905 (Fig. 4; Table S1 [suppl.]).

We note that natural disturbances are also important for regeneration of conifers in this forest, aside from human management history. For example, it is possible that the observed lower peak in the regeneration of Abies firma from 1946 to 1965 was caused by wind disturbance, which also peaked in the same period. When these trees regenerated, more than 20 years had passed since the last clear cut was conducted for broad-leaved trees in 1919. In the absence of canopy gaps, low light levels in the forest understory likely limited the regeneration of Abies firma . Kabaya (1975) reported saplings of Abies firma growing in natural gaps in a nearby fir-hemlock forest at UTCBF.

The frequency distributions of tree regeneration year for the two broad-leaf tree species ( Cinnamo­mum yabunikkei and Quercus acuta ) were consistent with the above-described management history, and also supported our hypothesis that regeneration was stimulated by wind-related disturbances. In this forest, understory broad-leaves were cut once every 20-30 years until the 1920s. Considering the time required for seedlings to reach a height of 1.37 m from germina­tion, the oldest regeneration year for broad-leaf trees observed in this forest (1923) corresponded well with the last recorded clear cut. The regeneration year of these species peaked from 1946 to 1965, consistent with a peak in the expected number of wind-damage events during this period (Figs. 4 and 5). After the end of coppice-with-standards management around 1920, regeneration of shade-tolerant broad-leaved species was likely dependent on canopy gaps formed by strong wind events. However, the temporal regeneration pattern of Castanopsis sieboldii could not be explained by management history or wind disturbance. Given that it is a shade-tolerant species, its regeneration may be less dependent on disturbance factors.

We developed a flexible model framework for the analysis of longitudinal wind-damage records for forests that adequately addresses autocorrelation. In an application example, we demonstrated the usefulness of our model for understanding forest regeneration history. The model presented in this study can be readily expanded to include other factors representing additional forest, wind, and environmental conditions. The models discussed here represent a simple framework for better understanding past forest dynamics and for making appropriate forest management plans to predict and mitigate wind damage.