The 30-year impact of post-windthrow management on the forest regeneration process in northern Japan

The frequency and intensity of typhoons are expected to increase over time due to climate change. These changes may expose forests to more windthrow in the future, and increasing the resilience of hemiboreal forests through forest management after windthrow is important. Here, we quantified forest structure recovery using aerial photos and light detection and ranging (LiDAR) data after catastrophic windthrow events. Our aims are to test the following three hypotheses: (1) forest structure will not recover within 30 years after windthrow, (2) forest recovery will be affected not only by salvaging but also pre-windthrow attributes and geographical features, and (3) various post-windthrow management including salvaging will drastically alter tree species composition and delay forest recovery. Our results revealed that hypothesis (1) and (2) were supported and (3) was partially supported. The ordination results suggested that more than 30 years were needed to recover canopy tree height after windthrow in hemiboreal forests in Hokkaido, Japan. Salvage logging did not delay natural succession, but it significantly decreased the cover ratio of conifer species sites (0.107 ± 0.023) compared with natural succession sites (0.310 ± 0.091). The higher the elevation, the steeper the site, and the higher the average canopy height before windthrow, the slower the recovery of forest stands after windthrow and salvaging. Scarification and planting after salvage logging significantly increased the number of canopy trees, but those sites differed completely in species composition from the old growth forests. Our study thus determined that the choice and intensity of post-disturbance management in hemiboreal forests should be carefully considered based on the management purpose and local characteristics.


Introduction
The severity and frequency of wind disturbances have intensified due to climate change (Usbeck et al. 2010;Gregow et al. 2017;Laapas et al. 2019). This intensification may result in existing forest ecosystems being transformed into different or non-forest states with lower capacities to provide desired ecosystem services in the future (Kamimura et al. 2008;Johnstone et al. 2016). Nayak and Takemi (2019) predicted an increased number of typhoons with more precipitation in Hokkaido, northern Japan, in the future, and Murakami et al. (2012) also projected forests will experience more intense winds in northern Japan. Furthermore, hemiboreal forests dominated by conifer species are reported to be more sensitive to windthrow than those dominated by broadleaf species (Rich et al. 2007).
With the aim of maintaining the ecosystem services provided by hemiboreal forests and increasing their resilience to different disturbances, monitoring forest structure recovery after windthrow is of greater importance than ever before. A complex forest structure (i.e., more tree size diversity and the development of multilayered forests) maintains biodiversity (Maltamo et al. 2005). Forest structure is also the dominant driver of timber productivity (Bohn and Huth 2017). Thus, monitoring forest structure recovery with the aim of maintaining the ecosystem services provided by hemiboreal forests is of paramount importance. However, windthrow damage to forest structures is easily ignored. Although forest structures after windthrow have been studied intensively in temperate forests, the long-term impacts of windthrow on boreal hemiboreal forests remain uncertain (Rich et al. 2007;Taeroe et al. 2019). Senf et al. (2019) found that forest structures on 84% of the disturbed area on average reached the recovery threshold within 30 years post-disturbance in temperate forests in central Europe. As the growth of trees in hemiboreal forests is relatively slower than that in temperate forests, the forest structure recovery after a windthrow event may also be slower. Kosugi et al. (2016) showed that more than 60 years are needed for hemiboreal forests to reach the stem exclusion stage following disturbance.
However, forest structure recovery can be altered by various factors. Post-windthrow management is a key factor affecting forest structure recovery after windthrow, and the potential damage associated with windthrow has been largely overlooked until very recently (Seidl et al. 2014;Taeroe et al. 2019). Salvage logging aims to prevent subsequent fires and insect outbreaks and can remove nurse logs and might destroy residual vegetation, the main sources of subsequent recovery (Morimoto et al. 2011). The consequences of these actions might last for 60 years (Morimoto et al. 2019). These practices also reduce the carbon stocks following catastrophic windthrow events (Hotta et al. 2020), potentially aggravating climate change in return. Soil scarification aims to remove organic-rich soil layers, which contain pathogens, and understory vegetation, which inhibits tree regeneration. However, advanced tree seedlings are also completely destroyed; thus, forest recovery would be much slower in scarified sites (Nilsson et al. 2006;Aoyama et al. 2011). Planting conifers is a traditional way to provide industrial wood, and conifers are more sensitive to windthrow. This process is thus likely to create forests that are more vulnerable to windthrow (Dhubháin and Farrelly 2018). In the face of future climate change, adaptation of the current management practices is urgently needed.
However, the impacts of forest management are mostly geographically confined. More attention should thus be given to local characteristics, such as the pre-windthrow attributes and geographic features that affect forest recovery (Barij et al. 2007;Hautier et al. 2018;Hu et al. 2018). The regeneration process following windthrow is based on attributes such as dispersed seeds and existing saplings from the pre-disturbance period in most conifer forests in the hemiboreal zone (Greene et al. 2011;Hérault and Piponiot 2018). Differences in pre-windthrow attributes could thus possibly affect successional composition (Kosugi et al. 2016). Slope aspects connected with light competition have also been suggested to be crucial in the early stage of forest recovery after a disturbance (Hasler 1982;Hu et al. 2018). Elevation can affect the structural complexity of a forest (Asner et al. 2014;Tudoroiu et al. 2016), and steepness might influence the forest structure by altering the water uptake and seed erosion processes (Barij et al. 2007). Considering these effects should be necessary when formulating suitable management strategies targeted at various regions.
Photogrammetry and LiDAR data have been proven to be powerful methods for assessing forest resources over a wide range (Shi et al. 2020). However, no previous studies have managed to use it to conduct landscape-level research, which obtains a combined understanding based on various types of major management practices (e.g., salvage logging, scarification, planting and seeding), predisturbance attributes and topographic features after windthrow (Rammig et al. 2007).
Here, we attempted to quantify the forest structure recovery 30 years following a major wind disturbance event using aerial photos and LiDAR data. We specifically examined the impacts of five types of postmanagement practices, different predisturbance attributes and topographic features on the forest structure recovery process in a hemiboreal forest area in Hokkaido, northern Japan, which plays an irreplaceable role in alleviating global warming. As our results may be utilized to optimize post-management processes in the face of increasing wind disturbances at the landscape level, we raised three hypotheses: (1) forest structure will not recover within 30 years after windthrow, (2) forest recovery will be affected not only by salvaging but also pre-windthrow attributes and geographical features, and (3) various postwindthrow management including salvaging will drastically alter tree species composition and delay forest recovery.

Study area
We targeted the University of Tokyo Hokkaido Forest (UTHF) in Furano, Hokkaido, Japan ( Fig. 1a; 43° 10′-21′ N, 142° 23′-41′ E). Various forest management experiments have been conducted and recorded for this area since 1958. The available aerial photos before and after windthrow and recent LiDAR data that can be found for that area enable this research. The total forest area is 22,717 ha, with a large elevational difference from 190 to 1459 m. From 2011 to 2020, the average temperature at the arboretum (230 m) was 6.6 °C, and the annual precipitation was nearly 1196 mm (The University of Tokyo Hokkaido Forest 2021). From the end of November to the beginning of April, the forest is covered by snow, and the maximum snow depth is 85.6 cm. Abies sachalinensis (F. Schmidt) Mast (Sakhalin fir) dominantly covers the UTHF site from elevations of 200 m to approximately 1200 (Jayathunga et al. 2018), and dwarf bamboo (Sasa senanensis (Franch. Et Sav.) Rehder and S. kurilensis (Rupr.) Makino et Shibata occupy the forest floor (Owari et al. 2011). A strong typhoon (Typhoon Thad) affected the UTHF on the 23rd of August, 1981. The 1 3 wind reached 26.3 ms −1 and caused severe windthrow within most of the forest area (Fig. 1b). After this typhoon, a recovery plan including 4 forms of post-windthrow management (Table 1)  To reveal the combined influences of the post-windthrow management projects, pre-windthrow attributes and topographic features on forest structure recovery, we selected 48 vegetation plots ( Fig. 1c 10 X X X -X (Oak species) SL_P(As) 10 X X X X (Abies species) -SL_P(Pg) 10 X X X X (Picea species) -1 3 50 m × 50 m based on a wind disturbance map (Fig. 1b) created in ArcGIS (ArcMap 10.7.1, Esri Japan Corp., Japan) by referring to management records from 1981 to 1986. All plots were focused on the seriously damaged area determined by assessing aerial photos taken in 1981. We also certified that no other major artificial or natural disturbances occurred in the study plots after the post-management project by checking the records from 1986 to 2015 provided by UTHF. The reference (RF) in our study was represented by a stand of forests that has undergone no major disturbances in the last 30 years. We considered this reference plot to represent the possible climax stage of the UTHF plots following a major wind disturbance event on a long-term scale.

Remote sensing data
In this study, we used aerial photos taken in 1977 to quantify the complexity of the canopy layers before the windthrow event and aerial photos taken in 2017 to quantify the complexity of the recent canopy layers 30 years following windthrow. Light detection and ranging (LiDAR) data (Table S1.1) collected from 2015 to 2019 were used to generate the leaf area index (LAI) data to determine the complexity of the current forest structure, including the shrub layers. Although a temporal difference was found between the acquisitions of the aerial photos taken in 2017 and the LiDAR data, we supposed that any influence on the results was negligible because (1) no major disturbance occurred in the study sites from 2015 to 2020 and (2) a 5-year gap is insufficient to produce any remarkable changes in hemiboreal forests. We collected aerial photographs (Table S1.2) targeting the study sites from the Ministry of Agriculture, Forestry and Fisheries in Japan (MAFF) and the Geospatial Information Authority of Japan (GIA). We quantified the average canopy tree height in 1977 as the pre-windthrow attributes and the average canopy tree height, number of canopy trees and forest cover ratio in 2017 as the forest structure 30 years after the management projects using Stereo Viewer Pro (PHOTEC Co., Ltd., Japan). Stereo Viewer Pro allows the user to perform stereoscopic measurements based on aerial photographs. It is an established method for measuring canopy tree height and forest cover ratio, and the validity of data using this method has been demonstrated in several studies (St-Onge et al. 2004;Korpela 2004;Korpela and Tokola 2006;Shimizu et al. 2014). In Stereo Viewer Pro, we measured the height of all visible vegetation in each plot ( Fig. S1.1). Each measured vegetation was set as a 3D point with X, Y, and Z coordinates in Stereo Viewer Pro; these points could then be opened in ArcGIS or exported as a text file. The polygon of each crown was created in Stereo Viewer Pro and transferred into ArcGIS for the following analyses. In ArcGIS, we first merged the crown polygons to represent the forested area in each plot and used the field calculator function to calculate the average canopy tree height, number of canopy trees and forest cover ratio in each plot (Table 2). We confirmed the accuracy of tree height estimated by Stereo Viewer Pro by comparing with tree height estimated by LiDAR ( Fig. S1.2). We found a strong consistency between them. Further, we compared the estimated canopy tree density with stem density derived from field census data, provided by UTHF, for 11 plots ( Fig. S1.3). The stem density from remote sensing was roughly consistent with the density of stems with DBH ≥ 11 cm from the field data. East, south, west and north. The east-facing slope included the slopes facing from northeast to southeast; the south-facing slope included the slopes facing from southeast to southwest; the west-facing slope included the slopes facing from southwest to northwest; and the northfacing slope included the slopes facing from northwest to northeast. Furthermore, we aimed to compensate for the limitation whereby aerial photos might neglect trees in the intermediate or suppressed layers by estimating the LAI using LiDAR data. The LiDAR datasets used herein were acquired in 2015, 2017, 2018 and 2019 using an Optech Airborne Laser Terrain Mapper (ALTM) Orion M300 sensor (Teledyne Technologies, Waterloo, ON, Canada) and a GNSS Trimble 5700・R3 (Trimble Inc., USA) mounted on an AS350B helicopter (Eurocopter Group, USA) (Moe et al. 2020). Table 4 shows the flight parameters of the LiDAR campaigns each year. The LiDAR data were initially processed and delivered in LAS format by Hokkaido Aero Asahi Corp. (Japan). The accuracy of these LiDAR data was also verified by Jayathunga et al. (2018) and Moe et al. (2020). The ALS Forest function in LiDAR360 V4.1 (GreenValley International Corp., USA) enabled us to calculate the LAI in each plot from the LiDAR data provided by the UTHF. In this process, we first normalized the point cloud data by using the 'Normalize by DEM' function in LiDAR360. Then, we used the 'ALS Forest' function to calculate the LAI by testing several parameters, including the grid size, height break and leaf angle distribution. We finally chose the recommended value of 0.5 as the leaf angle distribution (Li et al. 2016;Jin et al. 2020) and 20 m as the grid size. As the maximum individual crown in our field measured by Stereo Viewer Pro was larger than 15 and less than 20 m, we chose 20 m instead of the recommended 15 m when running the function. The height break was a parameter that allowed us to avoid the influence of vegetation on the forest floor. The final value we chose for this parameter was 3 m, as we considered any vegetation lower than 3 m to be representative of the shrub layer when processing the aerial photos previously.
We performed a field survey on the 3rd and 4th of August 2021 to identify the dominant species in the canopy layer. We first visually separated the identified trees into 237 groups using aerial photos based on texture and color and labeled these groups manually. However, considering the time required for this visual identification process, we selected only 6 plots under each management scheme for conducting the species classification process. In this process, we selected one representative tree from each of the 237 groups and circled these trees in both the aerial photos and digital canopy height model (DCM; resolution of 0.1 m) provided by UTHF, thus collecting the crown images and the GNSS data. We brought the prepared data to the field to verify the species of the targeted trees using aerial photos and a handy GNSS unit (GPSMAP ® 62SCJ, Garmin Ltd.). Twenty-one species represented by 237 trees were recorded in the field (Table 3). The subsequent analyses were performed at the species level basically. But four species of Acer, two species of Betula, and two species of Salix were analyzed at the genus level because they sometimes appear to be the same in aerial photographs. We also classified the species into three regeneration types based on Hanada et al. (2006): pioneer species, gap species, and shade-tolerant species.

Statistical analyses
R statistical software (v. 4.0.3, R Core Team 2020) was used for the statistical analysis conducted in this study.

Forest structure characteristics 30 years after the windthrow event and post-windthrow management practices: principal component analysis (PCA)
To determine the forest structure characteristics following the windthrow event and according to each post-windthrow management scheme, a PCA was performed. The standardized indicator data (the average canopy tree height, number of canopy trees, forest cover ratio and LAI data) were analyzed by PCA using the 'prcomp' function in R. The length of each arrow indicates the variance of the indicators in the ordinations, whereas the angles between the arrows indicate the strength of the correlation between the indicators (Gniazdowski 2017). The 95% confidence ellipse of each group indicates the overall characteristics of and similarity among groups. In the PCA results, if the 95% confidence ellipse was closer to the reference ellipse with a similar angle of inclination, we considered this sample to have a more similar structure with the reference, which might imply a more complex structure (Franklin and Van Pelt 2004).

Species composition characteristics 30 years after windthrow and post-windthrow management: pairwise test
To investigate the differences in the canopy-layer species composition among the management categories, we conducted pairwise tests for multiple comparisons. We ran the 'pairwise.t.test' function in R and chose Bonferroni as the method for adjusting the p values to compare the relative abundance of classified species in the canopy layer, including broadleaf, conifer, pioneer, and shade-tolerant species.

Domination of species and species composition similarity 30 years after the windthrow event and post-windthrow management: nonmetric multidimensional scaling (NMDS)
To diagrammatically determine the species composition similarity among the management groups, the species 1 3 a Gap species: a species whose shade tolerance is between those of pioneer species and shade-tolerant species composition was ordinated by NMDS using the 'meta-MDS' function in the package 'vegan' (version 2.5-7). In the NMDS process, we used the log-transformed coverage data of each post-windthrow management group and selected the Bray-Curtis scale as the length index.

Factors influencing forest structure recovery 30 years after windthrow: generalized linear models (GLMs)
To recognize the important factors affecting forest structure recovery at the landscape level, we designed 2 generalized linear models (GLMs) using the 'stats' package in R. To identify the determinants of the average canopy tree height in 2017, we chose the average canopy tree height in 2017 as the response variable and the post-windthrow management schemes (NS and SL), the average canopy tree height in 1977 (representing the pre-windthrow attributes) and topographical features as the explanatory variables in Model 1 (Table 4, Model 1; average tree height derived using the gamma distribution error with the inverse link function). We selected the average canopy tree height in 1977 in Model 1 as a factor because taller trees exhibit a higher windthrow risk (Krejci et al. 2018). A slow recovery will be easily found in serious windthrow areas. Thus, we expected that the higher average canopy tree height before windthrow, related to a more serious windthrow, would negatively affect the average canopy tree height recovery 30 years after windthrow.
To recognize the factors influencing forest structure recovery in 2017 following scarification, sowing and planting, we chose the average canopy tree height, number of canopy trees, forest cover ratio and LAI in 2017 as the response variables. We selected the post-windthrow management schemes (NS, SL, SL_S, SL_P(As) and SL_P(Pg)) and topographical features as the explanatory variables in Model 2 (Table 4, Model 2; average canopy tree height derived using the gamma distribution error with the inverse link function; number of canopy trees derived using the gamma distribution error with log link function; forest cover ratio derived using the Gaussian distribution after the arcsin transformation (arcsin(x^2)); and LAI derived using the gamma distribution error with the inverse link function). We eliminated the average canopy tree height in 1977 from the explanatory variables when choosing the post-windthrow management schemes (NS, SL, SL_S, SL_P(As) and SL_P(Pg)) as one of the explanatory variables because stand development was suggested by Hotta et al. (2020) to restart from "bare land" following scarification, thus thoroughly destroying advanced seedlings. We thus considered the limited influence caused by the average canopy tree height in 1977 after scarification. After the full model was built, candidate models were generated and ranked based on the ΔAICc values through the 'dredge' function in the package 'MuMIn' (version 1.43.47). All candidate models up to ΔAICc < 4 were extracted to calculate the average model using the 'model.avg' function following Anderson et al. (2001). We displayed the estimate, confidence interval (CI) (95% confidence limit) and relative value importance (RVI) of each explanatory variable from the full model-averaged coefficients derived in the results section. We concluded that a significant influence occurred when the CI of a continuous variable (i.e., the elevation or slope angle) excluded zero. We also used the 'emmean' function in the 'emmeans' package through Tukey's methods to compute the estimated marginal means (EMMs) to provide post hoc contrasts after model averaging. We used this approach to identify whether any significant differences existed among the categorical variables (i.e., the post-windthrow management scheme or slope aspect), as suggested by Sumasgutner et al. (2016).

Forest structures at various post-windthrow management sites
The PCA results (Fig. 2) revealed the forest structure of each study plot 30 years after the windthrow event. Only the first two axes were presented and interpreted due to these axes corresponding to a meaningful proportion (91.98%) of the explained variance. The first principal component (Comp. 1) was correlated with all indexes (average canopy tree height, number of canopy trees, forest cover ratio and LAI) informing the forest structure, whereas the second principal component (Comp. 2) was also correlated with all of these indexes except average canopy tree height. All reference forest (RF) plots were intensively plotted at the rightmost ordination and at the top of the y-axis. We considered higher average canopy tree height, forest cover ratio and LAI values and a lower number of canopy tree values to represent the RF. The windthrow sites mostly plot to the left of the RF. This implies that the windthrow sites correspond to lower canopy tree height and LAI values compared with the RF sites. The 95% confidence ellipse of NS was the nearest ellipse to that of RF. This result means that NS had the forest structure most similar to that of RF among all of the groups. On the other hand, SL_P(Pg) was the farthest from RF, showing a significantly different forest structure from that of RF.

Species composition of the canopy layers at various post-windthrow management sites
The species composition at various post-windthrow management sites among the 6 analyzed groups is shown in Fig.  S2.1. Both conifer species and broadleaf species covered nearly 50% of the canopy layer in the RF sites (Fig. 3a, b). Compared with the RF site, significant increases were found in conifer species in SL_P(As) and SL_P(Pg), while dramatic decreases in conifer species were found in SL and SL_S. No significant difference was observed between RF and NS in terms of the conifer or broadleaf cover ratio.
On the other hand, shade-tolerant species dominated the canopy layer of RF, which was related to the low pioneer species cover ratio (Fig. 3c, d). Among RF, NS, SL_P(As) and SL_P(Pg), the coverage ratios of pioneer species were almost the same, while the coverage ratio of shade-tolerant species was slightly lower at the NS site than at the other three sites. At the SL site, the number of pioneer species was higher, while that of shade-tolerant species was significantly lower than the corresponding values in the R and NS plots. In contrast, shade-tolerant species covered only 25% of the canopy layer, and the coverage ratio of pioneer species significantly increased in SL_S compared with the other sites. The NMDS results (Fig. 4) show the canopy-layer species composition at different post-windthrow management groups. As the stress value was 0.120, which was lower than 0.200, we assumed two dimensions would be acceptable when running the NMDS. The plots of the RF sites were located to the left of the x-axis, and the NS sites were adjacent to the x-axis. Although a significant difference in species composition was shown between the two groups, NS was more similar to RF than SL, SL_S or SL_P(Pg) ( Table S2.1). The results also explained that shade-tolerant conifer species, such as P. jezoensis and A. sachalinensis, tended to occur more in RF (Fig. S2.1). SL_S was intensively plotted above the y-axis at the same location as the pioneer species Betula. SL occurred to the bottom left and was located between NS and SL_S. Although the species composition among NS, SL and SL_S varied from group to group, SL was more similar to NS and SL_P(As) but differed from SL_S and RF (Table S2.1). The coverage of broadleaf species, especially T. japonica, was higher in SL (Figs. 4, S2.1). SL_P (As) plots were separately plotted at the bottom left of the ordinate. A. sachalinensis appeared to occur in this area, as shown by the NMDS results. SL_P(Pg) was the only group located to the right of the x-axis, and P. glehnii tended to occur more frequently at this site.

Factors affecting the regeneration process after windthrow
Analyzing the canopy tree height derived in 2017 demonstrated that slope aspect significantly affected canopy tree height recovery (Table 5). The predicted average canopy tree heights on the south side of the study plot were significantly higher than those on the east side, as revealed by the post hoc test (Fig. 5, Table S4.1). We also found that forests with higher average canopy tree heights before the windthrow event had lower average canopy tree heights 30 years after the windthrow event (Table 5). However, as the CIs of the post-windthrow management and the slope angle included 0, our results cannot support any significant influence of the analyzed post-windthrow management schemes (NS and SL) or the slope angle on the average canopy tree height recovery after the windthrow event in the natural succession or salvaged sites.
The results of Model 2 suggested that post-windthrow management had a strong influence on the regeneration process related to the average canopy tree height ( Table 6). The predicted values derived from the GLM and post hoc processes by estimating the marginal means (Fig. 6) showed no significant difference between NS and SL, while surprisingly higher stem densities were estimated by the GLM (Fig. 6b, Table S6.1) in SL_S, SL_P(As) and SL_P(Pg) compared with those of NS and SL. However, no clear impact was found on the average canopy tree height or LAI due to NS, SL, SL_S or SL_P(As) according to the GLM results (Fig. 6a, d). Although the estimated forest coverage ratios were the same in NS, SL and SL_S, these ratios were significantly lower in SL_P(As) and SL_P(Pg) (Fig. 6c, Table S7.1). The LAI values predicted by the GLM also substantially decreased in SL_P (Pg) compared with those derived for the other sites (Fig. 6d, Table S8.1).
At the same time, the number of canopy trees, forest cover ratio and LAI were also related to the post-windthrow management scheme (including scarification, planting and seeding). Our results revealed that the average canopy tree height was more likely to decrease at high elevations (Table 6). Similarly, the forest coverage ratio and LAI in 2017 were negatively influenced by elevation. However, although elevation was selected in the final model corresponding to the number of canopy trees, no significant influence was found, as the 95% confidence CI ranged from − 1.11E−06 to 1.44E−06. We concluded that the slope angle did not affect the average canopy tree height, forest cover ratio or LAI for the same reason, but the results implied that higher stem densities could be found more easily in flat landscapes than in steep landscapes (Table 6).

Thirty-year impact of catastrophic windthrow events
Our study suggested that the forest cover ratio underwent a full recovery approximately 30 years after the windthrow event. However, the LAI and average canopy tree height (Fig. 2) values of the disturbed areas were still significantly lower than those of the reference area. From these results, we believe that the current forests are still in the recovery process during which they are developing stems (Oliver 1980) even 30 years after the windthrow event. In other research (Götmark and Kiffer 2014) conducted in a hemiboreal zone in Sweden, the average recorded canopy tree height was nearly the same as that reported in this work even 40 years after a windthrow event. Although no significant differences were found in the coverage ratio of the functional types, the nonsalvaged area tended to have a higher coverage ratio of pioneers and broadleaved trees but a lower coverage ratio of shade-tolerant and coniferous trees than the reference (Fig. 3). In addition, we found a divergence of species composition between the nonsalvaged area and reference (Fig. 4). Conifers such as Picea spp. and A. sachalinensis usually dominate undisturbed forests in Hokkaido, northern Japan, while broadleaved trees such as Betula spp. become dominant in windthrow areas due to their fast growth rate in the gap environments (Oliver 1980). Furthermore, conifers could be more severely damaged by wind disturbance because conifers are more susceptible to wind disturbance than broadleaves (Dhubháin and Farrelly 2018). Most conifers might be destroyed during the wind disturbance and recover slowly due to the gap environments created by windthrow, which cause the divergence of species composition between the nonsalvaged and reference areas.

Impact of salvage logging, pre-windthrow attributes and geographic features on forest recovery rate
The forest structure recovery was not evidently changed after 30 years of salvage logging (Fig. 6), which was similar to several previous studies (Royo et al. 2010;Morimoto et al. 2011). However, we found a significant decrease in conifer species but an increase in pioneer species at the salvaged sites compared with the nonsalvaged area (Fig. 3). Similar results were also reported by Orczewska et al. (2019), who observed that the number and coverage of ancient forest indicator species decreased following salvage logging. Salvage logging often severely destroys residual vegetation, including tree saplings (Morimoto et al. 2011), which can largely affect the subsequent recovery process. Indeed, Kurahashi et al. (1986) reported that approximately 25% of saplings were destroyed by salvage logging after windthrow at our study site. However, several geographic features and the canopy tree height before windthrow actually affect the forest regeneration rate. We found that the slope aspect affected the growth of individual trees in the naturally regenerating area ( Table 5). The south-facing aspect, which corresponds to increased light availability in the Northern Hemisphere, facilitates the growth of light-demanding species in disturbed areas (Vodde et al. 2010). The slope aspect can also affect the abundance of advanced regeneration in forest floors; similarly, Owari (2013) showed that A. sachaliensis saplings were abundant on southwestern slopes. Regeneration from saplings enables fast recovery of canopy height after wind disturbance, which could cause greater heights on south-facing slopes (Fig. 5). Limited resources and low temperatures at high elevations (Coomes and Allen 2007;Potapov et al. 2020) can be the main reasons for the slow regeneration process observed at high elevations, as seed erosion occurring on steep slopes (Barij et al. 2007) accounted for a lower number of canopy trees compared with flat areas (Table 6). Furthermore, we found that the canopy tree height in 1977 had a negative influence on the canopy tree height in 2017 (Table 5). This result indicates that the average post-windthrow canopy heights tended to be higher in the forests where the average canopy tree height was lower before the windthrow event. As the windthrow  Table 6 for the model-averaged coefficients). The black triangles show the mean predicted values derived through model averaging. The gray circles and boxplots suggest the observed value (south or east, n = 18). Significant differences (p < 0.05) are shown between the groups with different letters according to the estimated marginal means (EMMs) derived using Tukey's method (see Table S3.1 for the p value) risk increases with the growth of the canopy tree height (Oliver 1980;Dhubháin and Farrelly 2018), the forests in which canopy heights were lower have the advantage of high wind resistance. The stands where more surviving trees remain due to the lower windthrow risk will attain a rapid recovery of average canopy tree height.

The impact of scarification, planting and sowing on species composition and factors affecting forest recovery rate
Clear increases in the coverage and abundance of pioneer species, such as Betula spp., were found in the sites that Table 6 Coefficients estimated by the full average model built using Model 2 for each metric related to the forest structure The models corresponding to average canopy tree height and LAI were built through gamma error with the inverse link function, and their estimated results are given on the inverse scale. The model corresponding to the number of canopy trees was built through gamma error with the log link function, and its estimated results are given on the log scale. The model corresponding to the forest cover ratio was built through Gaussian error with the identity link function, and its estimated results are given on the normal scale CI, confidence interval; RVI, relative variable importance; SL, salvage logging; SL_S, salvage logging + scarification + Quercus crispula sowing; SL_P(As), salvage logging + scarification + Abies sachalinensis planting; SL_P(Pg), salvage logging + scarification + Picea glehnii planting  (Figs. 3, 4). This indicates that the establishment of Q. crispula after scarification was clearly unsuccessful at our study site. The capacity of oak to become established might be reduced because scarification creates an increased light intensity and dry conditions (Wetzel and Burgess 2001;Asada et al. 2017), while Betula spp. gain an advantage as light-demanding species. Although previous research mentioned that scarification contributed to the successful establishment of oak following direct seeding (Birkedal et al. 2010), seeding Q. crispula after scarification might be unsuitable at our study site. Planting directly can lead to an identically high number of canopy trees recovering (Fig. 6b), while it might also hamper the natural regeneration process (Thorn et al. 2017). Planting can result in a huge decrease in species diversity (Newbold et al. 2015). Species composition of the canopy layer was shifted in the P. glehnii planting plots, which were mostly dominated by P. glehnii (Fig. 4) and the A. sachalinensis planting plots were dominated by planted A. sachalinensis, whereas the species composition of SL_P(As) was more similar to the reference and natural succession sites (Table S2.1). The coverage of unplanted species and LAI in the A. sachalinensis planting area also tended to be slightly higher than those in the P. glehnii planting area (Figs. 6d, S2.1). Okamura et al. (1995) reported that broadleaved species mostly regenerated on piles of debris, where no trees were planted, in the P. glehnii planting area, whereas many  Table 8 for the model-averaged coefficients. a Predicted average canopy tree height for the post-windthrow management sites in 2017. b Predicted number of canopy trees for the postwindthrow management sites in 2017. c Predicted forest coverage ratio for the post-windthrow management sites in 2017. d Predicted LAI for the post-windthrow management sites in 2017. The black triangles show the mean predicted values derived through model averaging. The circles and boxplots suggest the observed values (NS, SL, SL_S, SL_P(As), and SL_P(Pg); n = 48). Significant differences (p < 0.05) between the groups are shown with different letters according to the estimated marginal means (EMMs) derived using Tukey's method (see Tables S4.2, S5.2, S6.2, and S7.2 for the p value). NS, natural succession sites after the windthrow event in 1981; SL, salvage logging sites after the windthrow event in 1981; SL_S, salvage logging + scarification + Quercus crispula sowing sites after the windthrow event in 1981; SL_P(As), salvage logging + scarification + Abies sachalinensis planting sites after the windthrow event in 1981; SL_P(Pg), salvage logging + scarification + Picea glehnii planting sites after the windthrow event in 1981 1 3 broadleaved species invaded among planted trees in the A. sachalinensis planting area, suggesting that P. glehnii plantations are exclusive against other species. Although structure indicators except the number of canopy trees attained similar levels to NS and SL in the A. sachalinensis planting plots, those in the P. glehnii planting plots were much lower than those in the other plots (Fig. 6a, c, d). This can be because of the slow growth rate of P. glehnii.
Among all geographic factors, elevation in response to a harsh environment negatively affected the recovery rate after scarification and planting after salvage logging (Table 6). However, unlike in the natural succession area, the slope aspect was not selected as a significant factor in Model 2 ( Table 6). The possible reason for the difference is that the scarification has thoroughly removed residual vegetation, including advanced regeneration. Such homogeneous conditions created by scarification might make the effects of the slope aspect unclear.

Conclusions
The landscape-level forest recovery process 30 years following a windthrow event was identified in this study using aerial photos and LiDAR data. Although the canopy coverage ratios of the forests in the disturbed area underwent full recoveries, the forest sites were still in the stem-exclusion stage and required more time to undergo full height and LAI recoveries. No influence of salvage logging was found on the forest structure recovery process. However, we noticed a trend in which the composition of broadleaved tree species increased following salvage logging. Scarification increased the number of canopy trees but changed the species composition due to the regeneration of Betula species at a high value. Seeding had little effect after scarification because the oak species sown were not successfully established. Planting reduced the chance for other species to establish, thus resulting in a different species composition in the forest canopy. On the other hand, our study verified the influence of geographical factors and pre-windthrow attributes on the forest regeneration process. The south-facing aspect was found to have a positive impact on natural regeneration in hemiboreal forests, but this influence was limited in the scarified and planted areas. The high elevation and steep area of the forest sites led to slow regeneration after the windthrow event and management practices. The higher the average canopy tree height was before the windthrow event, the lower the mean height the forest recovered to after the windthrow event.
Our findings provide implications for management schemes in hemiboreal forests facing windthrow. We suggest passive restoration, 'doing nothing' after windthrow, to maintain the number of conifer species in hemiboreal forests. It is highly recommended in the south aspect area and was found to have a positive impact on natural regeneration. In addition, we also recommend passive restoration at high elevations and steep landscapes because of the high cost, limited accessibility and counterproductivity toward the recovery rate. On the other hand, it is suitable to conduct planting or scarification in flat areas with low elevations to obtain a number of commercial trees, such as Betula or Abies species, in a short time after windthrow. As the regeneration of Betula after the scarifcation will greatly impair the impact of sowing, we do not recommend sowing activities after scarification in our study area.
Due to the limitation of the method itself, we could only crudely estimate the LAI through LiDAR data and obtain a rough image of species composition using aerial photos. Further efforts, such as ground surveys or applying machine learning technologies, are necessary to obtain more accurate data with a high efficiency from aerial photos and LiDAR. We found an impact of pre-windthrow attributes, such as canopy tree height, on the regeneration process, possibly because it may affect the initial stage of regeneration directly after windthrow. However, we cannot prove this due to missing data. A long-term monitoring test is thus crucial to have both directly before and after windthrow events to confirm this relationship and determine further potential risks of applying post-windthrow management practices to hemiboreal forests in Hokkaido, northern Japan.