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

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. https://doi.org/10.5424/ fs/2019283-15200 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. Competing interests: The authors have declared that no competing interests exist. Correspondence should be addressed to Kiyoshi Umeki: umeki@faclty.chiba-u.jp Forest Systems 28 (3), e019, 12 pages (2019) eISSN: 2171-9845 https://doi.org/10.5424/fs/2019283-15200 Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria (INIA) 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


Introduction
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; 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 winddamage 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 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 longitudinal 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.
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 December 2019 • Volume 28 • Issue 3 • e019 3 data obtained from tree rings from an old-growth firhemlock (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.

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 winddamage probability and influencing factors) for model 1 are as follows: where D i 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(p i ); p i 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 (ve i ; continuous variable expressed in m s -1 ), mean stand volume (vo i ; continuous variable expressed in m 3 ha -1 ) for a decade (de i ; discrete variable), wind direction (di i ; discrete variable), and season (s i ; discrete variable); e 1 -e 2 are functions describing the effects of wind direction and season, respectively, whose values are autocorrelated with their neighboring values (explained below in detail); b 0 is an intercept term; b 1 and b 2 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 phenology, 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 1 and/or e 2 would be autocorrelated with those of neighboring levels (e.g. e 1 (NE) would be correlated with e 1 (NNE) and e 1 (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 e il is e i (i = 1 or 2) for the lth level; δ l is a set comprising neighboring levels; and σ ei 2 is the conditional variance. e il follows a normal distribution N, with the expected value being the average of the values for the neighboring levels. Note that interdependency among e il is circular for wind direction and season; the first level is a neighbor of the last level and vice versa. Because e i (i = 1 or 2) is non-identifiable, we added a constraint that ∑ l e il = 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 decadespecific mean stand volume (vo(de i )) were estimated assuming that where vo obs is the observed mean stand volume for a decade de j for which the record of mean stand volume was available, and σ vo 2 is a variance term expressing the degree of accuracy of the records. We also assumed that vo(de i ) followed a one-dimensional random walk model: where σ vo,de 2 is a variance term expressing auto correlation strength between mean stand volumes. The likelihood function of model 1 consisted of Equations 1 and 4, and σ vo,de 2 were estimated simultaneously with other model parameters.
In model 2, wind-damage probability is a function of wind speed, wind direction, season, and decade: where e 3 is a function describing the effect of decades and following a one-dimensional random walk model: where σ e3 2 is a variance term expressing the autocorrelation strength. e 3 acts as an intercept term in equation 5.

Model application
To provide an example of applying the models described above, we reanalyzed longitudinal winddamage 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 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 (3) A Bayesian model relating wind-damage probability to biotic and abiotic factors

Forest Systems
December 2019 • Volume 28 • Issue 3 • e019 5 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 longterm 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 (di i ) into 16 levels (N, NNE, NE, ENE,… NNW), season (s i ) into 12 levels (closely corresponding to months), and decade (de i ) 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 calculated 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 winddamage 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 method (Honda, 1912(Honda, , 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 clearcut of broad-leaved trees was conducted in 1919. A detailed description of this forest is given in Abrams et al. (2017).
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 regeneration 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 b 2 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 3 (de i ); Fig.  S4c [suppl.]). The Bayesian credible interval for winddamage 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
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.]).

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

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 clarified relevant factors determining the probability of wind-damage.
We introduced temporally or directionally autocorrelated terms (e 1 , e 2 , and vo in eqn. 2 for model 1; e 1 -e 3 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 3 (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.

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 winddamage 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 estimated the temporal change in mean stand volume and detected a positive effect of mean stand volume on winddamage 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 3 ) 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 constant in a population, the frequency distribution of Kiyoshi Umeki, Marc D. Abrams, Keisuke  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 firhemlock 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 shadetolerant 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; 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 broadleaved 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 (Cinnamomum 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 germination, 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.