A model for longitudinal data sets relating wind-damage probability to biotic and abiotic factors: a Bayesian approach

Kiyoshi Umeki (Umeki, K)

Graduate School of Horticulture, Chiba University, 648 Matsudo, Matsudo, Chiba, 271-8510 Japan.

Marc D. Abrams (Abrams, MD)

307 Forest Resources Building, Department of Ecosystem Science and Management, Penn State University, University Park, PA 16802 USA.

Keisuke Toyama (Toyama, K)

The University of Tokyo Chiba Forest, 770 Amatsu, Kamogawa, Chiba, 299‒5503 Japan.

Eri Nabeshima (Nabeshima, E)

Faculty of Agriculture, Ehime University, Tarumi, Matsuyama, Ehime, 790-8566 Japan.



Aim of study: To develop a statistical model framework to analyze longitudinal wind-damage records while accounting for autocorrelation, and to demonstrate the usefulness of the model in understanding the regeneration process of a natural forest.

Area of study: University of Tokyo Chiba Forest (UTCBF), southern Boso peninsula, Japan.

Material and methods: We used the proposed model framework with wind-damage records from UTCBF and wind metrics (speed, direction, season, and mean stand volume) from 1905–1985 to develop a model predicting wind-damage probability for the study area. Using the resultant model, we calculated past wind-damage probabilities for UTCBF. We then compared these past probabilities with the regeneration history of major species, estimated from ring records, in an old-growth fir–hemlock forest at UTCBF.

Main results: Wind-damage probability was influenced by wind speed, direction, and mean stand volume. The temporal pattern in the expected number of wind-damage events was similar to that of evergreen broad-leaf regeneration in the old-growth fir–hemlock forest, indicating that these species regenerated after major wind disturbances.

Research highlights: The model framework presented in this study can accommodate data with temporal interdependencies, and the resultant model can predict past and future patterns in wind disturbances. Thus, we have provided a basic model framework that allows for better understanding of past forest dynamics and appropriate future management planning.

Additional Keywords: dendrochronology; tree regeneration; wind-damage probability model; wind disturbance.

Abbreviations used: intrinsic CAR model (intrinsic conditional autoregressive model); MCMC (Markov chain Monte Carlo); 16 compass points = N, NNE, NE, ENE, E, ESE, SE, SSE, S, SSW, SW, WSW, W, WNW, NW, NNW (north, north-northeast, northeast, east-northeast, east, east-southeast, southeast, south-southeast, south, south-southwest, southwest, west-southwest, west, west-northwest, northwest, north-northwest, respectively); UTCBF (the University of Tokyo Chiba Forest).

Authors' contributions: Conceived and designed the study: KU and MDA. Collected data: KU, MDA, KT, and EN. Analyzed data: KU. Wrote the paper: KU. Revised the manuscript: MDA, KT, and EN.

Citation: Umeki, K., Abrams, M.D., Toyama, K., Nabeshima, E. (2019). A model for longitudinal data sets relating wind-damage probability to biotic and abiotic factors: a Bayesian approach. Forest Systems, Volume 28, Issue 3, e019.

Supplementary material: (Tables S1, S2; Figures S1–S5 and Appendices A1, A2) accompany the paper on FS’s website.

Received: 23 May 2019. Accepted: 15 Nov 2019.

Copyright © 2019 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International (CC-by 4.0) License.

Funding agencies/Institutions Project / Grant
Ministry of Education, Culture, Sports, Science, and Technology of Japan 23380085, 23380079, and 26450205
Japan Society for the Promotion of Science (to the second author) Guest Researcher Award

Competing interests: The authors have declared that no competing interests exist.

Correspondence should be addressed to Kiyoshi Umeki:





Material and Methods






Wind is a major disturbance in the forests of the Japanese archipelago and in other parts of the world (e.g. Nakashizuka & Yamamoto, 1987; Peterson, 2000; Svoboda et al., 2014). Wind causes damage to established trees and provides opportunities for trees to regenerate in natural forests (Abrams & Orwig, 1996; Abrams et al., 1999). Many important features of natural forests, including structure, diversity, and dynamics, are influenced by wind disturbances (e.g. Ulanova, 2000; Mitchell, 2013). However, these disturbances cause significant economic losses for forest managers (e.g. Peterson, 2000). Therefore, understanding wind disturbance in a given forest is key to ecologically and economically sustainable management. Furthermore, the frequency and severity of catastrophic wind events (typhoons, hurricanes, tornados, derecho, etc.) have increased and are predicted to increase further as a result of climate change (Banholzer et al., 2014; Gregow et al. 2017). Thus, it is critical that scientists and land managers have a comprehensive understanding of how significant wind events effect forests in various global regions.

The first step toward mitigating wind damage is understanding the relationships between wind-damage probability and influential biotic and abiotic factors (e.g. Gardiner et al., 2008; Kamimura et al., 2008; Hanewinkel et al., 2008; Albrecht et al. 2012). Previous studies have developed models that predict wind-damage risk and can be used to minimize or mitigate serious wind damage. These models can be divided into two classes, each with its own advantages and disadvantages: mechanistic models (e.g. Peltola et al., 1999, Gardiner et al., 2000, Gardiner et al., 2008) and statistical models (e.g. Hanewinkel et al., 2008, Albrecht et al., 2012).

Statistical models, the class to which the model presented in this study belongs, are easier to develop and implement than are mechanistic models. They generally require less data and information on the mechanisms of wind damage than do mechanistic models (i.e. mechanistic models require many input parameters). However, low applicability and, at times, low prediction accuracy are disadvantages of statistical models. These disadvantages can be caused by the use of inappropriate models, resulting in misleading conclusions and inaccurate prediction models (i.e. the statistical model combined with the estimated parameters) even when data are sufficient. Secondly, the use of limited data can lead to a lack of applicability over space and time for statistical models. For example, a model developed from data obtained from a single wind event may be a poor predictor of damage caused by another wind event (Kamimura et al., 2016).

One of the objectives of this study was to address the first potential cause of inadequate statistical models, i.e. inappropriate models. We present a statistical model that can analyze longitudinal wind-damage data obtained from fixed locations, thereby appropriately addressing autocorrelation. Autocorrelation can result in an inflated sample size and underestimated variance, exaggerating the effects of explanatory variables (Lichstein et al., 2002). Wind-damage data have been analyzed with various statistical models, including artificial neural networks (Hanewinkel, 2005), linear regression (Schütz et al., 2006), the weights of evidence method (Schindler et al., 2009), logistic regression (Schindler et al., 2009; Valinger & Fridman, 2011), generalized linear mixed models (Hanewinkel et al., 2008; Albrecht et al., 2012; Donis et al., 2018), generalized additive models (Schmidt et al., 2010), and classification and regression trees (Albrecht et al., 2012). However, autocorrelation was not appropriately addressed in most cases. For example, although autocorrelation in the residuals was assessed, it was not considered in parameter estimation (Hanewinkel et al., 2008; Schindler et al., 2009). Hanewinkel et al. (2008) found considerable temporal autocorrelation in the residuals of logistic regression models, modeled the autocorrelation with autoregressive techniques, and added a resultant autocorrelation component to their prediction models, but autocorrelation was not considered in parameter estimation in the logistic models. Schindler et al. (2009) assessed spatial autocorrelation in the residuals of their logistic model, but did not modify their prediction model because the spatial autocorrelation was deemed minor. A few studies have introduced terms to detect spatial autocorrelation in damage severity (Martin-Alcon et al., 2010; Hanewinkel et al., 2014).

The model presented here shares the second disadvantage associated with many statistical models, namely limited data. If the data used to develop a model are not comprehensive, the resulting predictions are generally inaccurate for conditions that differ from those under which the data were obtained. In our case, we analyzed longitudinal data from fixed locations, so their predictive accuracy is likely to be low when applied to other locations. However, accuracy is not problematic when this model is used to make predictions across time for the target locations. Moreover, this model may partially overcome the hindrance of limited data, as it can estimate inaccurate or incomplete data while estimating other model parameters.

The primary purpose of this study was to present a flexible model framework for analyzing longitu­dinal wind-damage data that adequately addresses autocorrelation and inaccurate or incomplete data. This framework is particularly useful for long-term forest datasets, such as those maintained by universities, institutes, and federal or local governments (Nilsson et al., 2004; Usbeck et al., 2010; Brázdil et al., 2018). In an example of model application, we reanalyzed an 85-year wind-damage record from the University of Tokyo Chiba Forest (UTCBF; Fig. 1), southern Boso Peninsula, Japan (Abrams et al., 2017), and developed a model that predicted wind-damage probability for the entire forest area.

Figure 1. Location of the University of Tokyo Chiba Forest (UTCBF) and study plots. Points in the lower right diagram represent plots. The interval of contours was 10 m.

Our secondary objective was to demonstrate the usefulness of our model in understanding regeneration processes in a natural forest. We estimated the temporal pattern in the number of wind-damage events during an 80-year period using the model for UTCBF and compared these predicted estimates to tree regeneration data obtained from tree rings from an old-growth fir–hemlock (Abies-Tsuga) forest found at UTCBF. Abrams et al. (2017) described detailed spatial and temporal regeneration patterns in the dominant conifers in this forest stand. They concluded that the coniferous species were less shade tolerant than the codominant evergreen broad-leaved species and that conifer regeneration was dependent on human-caused, rather than wind-caused, disturbances (i.e. forest management). In contrast, evergreen broad-leaved species in this stand were more shade tolerant, and their seedlings could survive under the canopy of coniferous trees (Honda, 1927; Kaji, 1975; Abrams et al., 2017). However, temporal regeneration patterns and the relationships between regeneration and natural or anthropogenic disturbance are not clear for evergreen broad-leaved species.

Material and MethodsTop

Wind-damage probability model

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:

Di ~ Bern (pi) (1)

logit(pi) = b0 + b1vei + b2vo(dei) + e1(dii) + e2(si), (2)

where Di is a binary variable taking a value of 1 if wind damage occurred in the target stands during the ith wind event and 0 otherwise, and is assumed to follow a Bernoulli distribution Bern(pi); pi is the probability that the ith wind event will cause damage in the target stands; conditions during a wind event are described by wind speed (vei; continuous variable expressed in m s-1), mean stand volume (voi; continuous variable expressed in m3 ha-1) for a decade (dei; discrete variable), wind direction (dii; discrete variable), and season (si; discrete variable); e1e2 are functions describing the effects of wind direction and season, respectively, whose values are autocorrelated with their neighboring values (explained below in detail); b0 is an intercept term; b1 and b2 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 e1 and/or e2 would be autocorrelated with those of neighboring levels (e.g. e1(NE) would be correlated with e1(NNE) and e1(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:

where eil is ei (i = 1 or 2) for the lth level; δl is a set comprising neighboring levels; and σei2 is the con­ditional variance. eil follows a normal distribution N, with the expected value being the average of the values for the neighboring levels. Note that interdependency among eil is circular for wind direction and season; the first level is a neighbor of the last level and vice versa. Because ei (i = 1 or 2) is non-identifiable, we added a constraint that ∑l eil = 0. Intrinsic CAR models have often been used to represent spatially autocorrelated random effects, but their application to temporal or directional data has been rare in forestry and ecology (e.g. Shaddick & Wakefield, 2002, Mizusaki et al., 2015). A feature of the intrinsic CAR model is that data are not needed for all levels (e.g. months) to estimate the expected values for each level. However, complete datasets do improve model accuracy.

Records for mean stand volume from UTCBF were imperfect and inaccurate (records were available for seven of nine decades); therefore, values for decade-specific mean stand volume (vo(dei)) were estimated assuming that

voobs(dej) ~ N(vo(dej), σvo2), (4)

where voobs is the observed mean stand volume for a decade dej for which the record of mean stand volume was available, and σvo2 is a variance term expressing the degree of accuracy of the records. We also assumed that vo(dei) followed a one-dimensional random walk model:

vo(dei) ~ N(vo(dei-1), σvo,de2),

where σvo,de2 is a variance term expressing auto­cor­relation strength between mean stand volumes. The like­lihood function of model 1 consisted of Equations 1 and 4, and σvo,de2 were estimated simultaneously with other model parameters.

In model 2, wind-damage probability is a function of wind speed, wind direction, season, and decade:

logit(pi) = b1vei + e1(dii) + e2(si) + e3(dei), (5)

where e3 is a function describing the effect of decades and following a one-dimensional random walk model:

e3(dei) ~ N(e3(dei-1), σe32)

where σe32 is a variance term expressing the auto­correlation strength. e3 acts as an intercept term in equation 5.

Model application

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 (, 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.

Dendrochronological sampling

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


Wind-damage probability model

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 b2 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 (e3(dei); 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 observed wind-damage events and predicted wind-damage probability

Figure 2. Predicted wind-damage probability for UTCBF related to (a) wind speed, (b) wind directions, (c) mean stand volume, and (d) seasons. Solid lines indicate posterior means, and dark and light grey colors indicate 80% and 95% credible intervals, respectively. Predictions were made by model 1.

Figure 3. Temporal change in the observed and estimated mean stand volume for UTCBF. Points indicate the observed mean stand volume. Solid lines indicate posterior means of the estimated mean stand volume, and dark and light grey colors indicate 80% and 95% credible intervals, respectively.

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

Figure 4. Temporal change in the observed (a) and expected (b) number of days with wind-damage events in forests in UTCBF. Different shades of grey indicate the contribution of wind events with different levels of wind-damage probability. Expected number of days with wind-damage events was calculated with model 1.

Tree regeneration in a fir–hemlock forest at UTCBF

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.

Figure 5. Frequency distribution of the regeneration year of five major tree species in a fir-hemlock forest in UTCBF.


Wind-damage model

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 (e1, e2, and vo in eqn. 2 for model 1; e1e3 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. e3(dei) 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.

Wind-damage probability in UTCBF

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 (e3) 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.

Regeneration patterns in a fir–hemlock forest at UTCBF

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.


We thank Chris Bouma for measuring tree age on the tree cores obtained for this study. We thank UTCBF for giving us permission to sample tree cores. We also thank members of UTCBF and students of Chiba University for their help in fieldwork.


Abrams MD, Orwig DA, 1996. A 300-year history of disturbance and canopy recruitment for co-occurring white pine and hemlock on the Allegheny Plateau, U.S.A. J Ecol 84: 353-363.

Abrams MD, Umeki K, Bouma C, Nabeshima E, Toyama K, 2017. A dendroecological analysis of past dynamics for a fir-hemlock forest on the Boso peninsula, southeastern Japan. Tree-Ring Res 73: 59-74.

Abrams MD, Copenheaver CA, Terazawa K, Umeki K, Takiya M, Akashi N, 1999. A 370-year dendroecological history of an old-growth Quercus-Abies-Acer forest in northern Japan. Can J For Res 29: 1891-1899.

Albrecht A, Hanewinkel M, Bauhus J, Kohnle U, 2012. How does silviculture affect storm damage in forests of south-western Germany? Results from empirical modeling based on long-term observations. Eur J For Res 131: 229-47.

Banholzer S, Kossin J, Donner S, 2014. The impact of climate change on natural disasters. In: Reducing disaster: Early warning systems for climate change; Zommers Z, Singh A (eds) . pp: 21-49. Springer, The Netherlands.

Besag J, York J, Mollié A, 1991. Bayesian image restoration with two applications in spatial statistics (with discussion). Ann Inst Statist Math 43: 1-59.

Brázdil R, Stucki P, Szabó P, Řezníčková L, Dolák L, Dobrovolný P, Tolasz R, Kotyza O, Chromá K, Suchánková S, 2018. Windstorms and forest disturbances in the Czech Lands: 1801-2015. Agr For Meteorol 250-251: 47-63.

Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, Brubaker M, Guo J, Li P, Riddell A, 2017. Stan: A probabilistic programming language. J Stat Softw 76: 1-32.

Choshi local meteorological observatory, 1969. Meteorological disasters in Chiba prefecture. Choshi local meteorological observatory, Choshi, Japam [in Japanese]

Choshi local meteorological observatory, 1987. Meteorological disasters in Chiba prefecture from 1969 to 1985. Japan weather association, Tokyo, Japam [in Japanese]

Donis J, Kitenberga M, Snepsts G, Dubrovskis E, Jansons A, 2018. Factors affecting windstorm damage at the stand level in hemiboreal forests in Latvia: case study of 2005 winter storm. Silva Fenn 52: article id 10009.

Ennos AR, 1997. Wind as an ecological factor. Trends Ecol Evol 12: 108-111.

Gardiner B, Peltola H, Kellomäki S, 2000. Comparison of two models for predicting the critical wind speeds required to damage coniferous trees. Ecol Model 129: 1-23.

Gardiner B, Byrne K, Hale S, Kamimura K, Mitchell SJ, Peltola H, Ruel J-C, 2008. A review of mechanistic modelling of wind damage risk to forests. Forestry 81: 447-463.

Gregow H, Laaksonen A, Alper ME, 2017. Increasing large scale windstorm damage in Western, Central and Northern European forests, 1951-2010. Sci Rep-UK 7: 46397.

Hanewinkel M, 2005. Neural networks for assessing the risk of windthrow on the forest division level: a case study in southwest Germany. Eur J For Res 124: 243-249.

Hanewinkel M, Breidenbach J, Neeff T, Kublin E, 2008. Seventy-seven years of natural disturbances in a mountain forest area - the influence of storm, snow, and insect damage analysed with a long-term time series. Can J For Res 38: 2249-61.

Hanewinkel M, Kuhn T, Bugmann H, Lanz A, Brang P, 2014. Vulnerability of uneven-aged forests to storm damage. Forestry 87: 525-34.

Honda S, 1912. The theory of Japanese forest zones. Miura Shoten, Tokyo, Japan. [in Japanese]

Honda S, 1927. Guidance of actual practice of aforestation. Miura Shoten, Tokyo, Japan. [in Japanese]

Joseph M, 2016. Exact sparse CAR models in Stan.

Kabaya H, 1975. Ecological studies on the vegetation of the Boso Mountains I. Distribution and structure of the natural fir-hemlock forests. Bull Tokyo Univ For 67: 51-62. [in Japanese with English summary]

Kaji M, 1975. Studies on the ecological status of Abies firma forest in the Boso peninsula. Bulletin of the Tokyo University Forests 68: 1-23.

Kamimura K, Gardiner B, Kato A, Hiroshima T, Shiraishi N, 2008. Developing a decision support approach to reduce wind damage risk - a case study on sugi (Cryptomeria japonica (L.f.) D.Don) forests in Japan. Forestry 81: 429-445.

Kamimura K, Gardiner B, Dupont S, Guyon D, Meredieu C, 2016. Mechanistic and statistical approaches to predicting wind damage to individual maritime pine (Pinus pinaster) trees in forests. Can J Fort Res 46: 88-100.

Lichstein JW, Simons TR, Shriner SA, Franzreb KE, 2002. Spatial autocorrelation and autoregressive models in ecology. Ecological Monographs 72: 445-463.[0445:SAAAMI]2.0.CO;2

Martin-Alcon S, Gonzalez-Olabarria JR, Coll L, 2010. Wind and snow damage in the Pyrenees pine forests: effect of stand attributes and location. Silva Fennica, 44: 399-410.

Mitchell SJ, 2013. Wind as a natural disturbance agent in forests: a synthesis. Forestry 86: 147-157.

Mizusaki D, Umeki K, Honjo T, 2015. Disentangling long- and short-term changes in perennial organ functions in seasonal environments: A model of foliar chlorophyll and nitrogen in saplings of four evergreen broad-leaved trees. Photosynthetica 53: 356-368.

Nakashizuka T, Yamamoto S, 1987. Natural disturbance and stability of forest community. Jan J Ecol 37: 19-30. [in Japanese with English summary]

Negisi K, 1997. Chronological notes of Tokyo University Forest in Chiba (Supplement): Happening found in the communications of the Past 4. Tokyo Univ For 36: 138-229. [in Japanese]

Nilsson C, Stjernquist I, Bärring L, Schlyter P, Jönsson AM, Samuelsson H, 2004. Recorded storm damage in Swedish forests 1901-2000. For Ecol Manag 199: 165-173.

Peltola H, Kellomäki S, Väisänen H, Ikonen V-P, 1999. A mechanistic model for assessing the risk of wind and snow damage to single trees and stands of Scots pine, Norway spruce, and birch. Can J For Res 29: 647-661.

Peterson C, 2000. Catastrophic wind damage to North American forests and the potential impact of climate change. Sci Total Environ 262: 287-311.

R Core Team, 2019. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL

Research Section of Tokyo University Forests, Tokyo University Forest in Chiba, 1974. Chronological notes of Tokyo University Forest in Chiba (1) Tokyo Univ For 18: 9-28. [in Japanese]

Stan Development Team, 2018. RStan: the R interface to Stan. R package version 2.17.3.

Ruel J-C, Pin D, Cooper K, 1998. Effect of topography on wind behaviour in a complex terrain. Forestry 71: 261-265.

Schindler D, Grebhan K, Albrecht A, Schönborn J, 2009. Modelling the wind damage probability in forests in Southwestern Germany for the 1999 winter storm 'Lothar'. Int J Biometeorol 53: 543-54.

Schütz JP, Götz M, Schmid W, Mandallaz D, 2006. Vulnerability of spruce (Picea abies) and beech (Fagus sylvatica) forest stands to storms and consequences for silviculture. Eur J For Res 125: 291-302.

Schmidt M, Hanewinkel M, Kändler G, Kublin E, Kohnle U, 2010. An inventory-based approach for modeling single-tree storm damage - experiences with the winter storm of 1999 in southwestern Germany. Can J For Res 40: 1636-52.

Shaddick G, Wakefield J, 2002. Modelling daily multivariate pollutant data at multiple sites. J Roy Stat Soc C-App 51: 351-372.

Svoboda M, Janda P, Bače R, Fraver S, Nagel TA, Rejzek J, Mikoláš M, Douda J, Boublík K, Šamonil P, et al., 2014. Landscape-level variability in historical disturbance in primary Picea abies mountain forests of the Eastern Carpathians, Romania. J Veg Sci 25: 386-401.

Ulanova NG, 2000. The effects of windthrow on forests at different spatial scales: a review. For Ecol Manag 135: 155-167.

Usbeck T, Wohlgemuth T, Dobbertin M, Pfister C, Bürgi A, Rebetez M, 2010. Increasing storm damage to forests in Switzerland from 1858 to 2007. Agr For Meteorol 150: 47-55.

Valinger E, Fridman J, 2011. Factors affecting the probability of windthrow at stand level as a result of Gudrun winter storm in southern Sweden. For Ecol Manag 262: 398-403.

Watanabe S, 2010. Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. J. Mach. Learn. Res. 11: 3571-3594.