Study landscape
The study landscape is a 12,169 ha hemi-boreal forest landscape located in the University of Tokyo Hokkaido Forests (UTHF; 43°10′-21′N, 142°23′-41′E; 350-1,000 m above sea level (a.s.l); Fig. 1). The annual mean temperature and precipitation at the Rokugo meteorological observatory (43°18′6′′N, 142°31′18′′E, 315 m a.s.l.), which is the closest site to the study landscape, are 5.5℃ and 972.6 mm, respectively (Japanese Meteorological Agency, 2012, representing the averages from 1981-2010). The study landscape is covered by hemi-boreal mixed forests dominated by conifers such as Abies sachalinensis (F. Schmidt) Mast., Picea jezoensis (Siebold et Zucc.) Carriére var. jesoensis; deciduous broadleaves such as Tilia japonica (Miq.) Simonk. and Betula ermanii Cham.; and herbaceous species such as Sasa senanensis (Franch. et Sav.) Rehder. The dominant soil type is mostly Andosols (parent material: andesite, rhyolite, volcanic ash, or dacite) with some Cambisols (parent material: rhyolite or dacite) (IUSS Working Group WRB 2015). The dominant natural disturbance of the study landscape is stand-replacing windthrow. There are records of windthrow caused by Typhoon Marie in 1954 and Typhoon Thad in 1981 in UTHF (Watanabe et al. 1990), and 8,735 ha (38.9% of the UTHF) in the area and 807,000 m3 in timber volume were damaged by Typhoon Thad in 1981 (Takada et al. 1986).
Simulation model and its application in this study
In this study, we applied LANDIS-II (Scheller et al. 2007), a forest landscape model, to predict the forest dynamics associated with climate change, changes in windthrow regimes, and post-windthrow management at the landscape scale. The landscape in LANDIS-II was represented as a collection of uniformly sized grid cells that contained the data for vegetation, climate, and soil properties. In this study, the size of a grid cell was defined as 1 ha. The vegetation data were represented as species-age cohorts, and each grid cell could contain multiple cohorts. The vegetation and environments influenced those of the other surrounding grid cells through seed dispersal. LANDIS-II requires a single succession extension and can include multiple disturbance or output extensions.
The improved version of LANDIS-II Net Ecosystem Carbon and Nitrogen succession v6.3 (Scheller et al. 2011; NECN succession) with tree-grass competition (Scheller et al. 2021) and the regeneration process on downed logs (Hotta et al. 2021) was used as the succession extension to simulate the dynamics of forest species composition and carbon balance. In NECN succession, biomass growth and seedling establishment are calculated based on the environmental conditions in each grid cell. The biomass growth of each cohort was calculated as the difference between monthly aboveground net primary production (NPP) and monthly mortality. Monthly aboveground NPP is calculated by multiplying the maximum value of monthly aboveground NPP (user-defined parameter for each species) by coefficients related to environmental limitation factors, such as monthly mean temperature, plant-available soil nitrogen, soil water content, the leaf area index (LAI) of the cohort itself, and the LAI of the other cohorts within the grid cell. Regarding cohort establishment, the establishment probabilities of each species in each grid cell were calculated based on the following two categories: (1) minimum January temperature (the minimum temperature of the coldest month in a year), growing degree days, and soil moisture content of each region; and (2) light availability at the site. If these two requirements were satisfied, cohorts were established in the grid cell. Additionally, tree species requiring dead wood for regeneration could be established only in grid cells that had well-decayed downed logs (Hotta et al. 2021). Two species of Picea spp., including P. jezoensis and P. glehnii, were set as dead wood-dependent species in this study. The decomposition rates of dead wood and soil organic matter were controlled by the soil temperature and moisture content. The decomposition rate of dead wood was calculated according to the amount of dead wood of each tree species in the dead wood pool using the parameters of dead wood decomposition rates of each species (user-defined parameter; WoodDecayRate). The decomposition rate of soil organic matter was calculated based on the decomposition rate parameters of the four soil pools (user-defined). Both decomposition rates were calculated to reflect the environmental conditions of each grid cell through soil temperature and moisture contents (Scheller et al. 2011).
LANDIS-II has been applied to various landscapes and ecosystems in North America, East Asia, Europe, and so on (e.g., Haga et al. 2020; Hotta et al. 2021; Lucash et al. 2019; Shifley et al. 2017). The following calibrations and validations in this study landscape were already reported in Hotta et al. (2021): the calibrations of aboveground biomass growth and litterfall of each species, the LAI and NPP in a grid cell, the aboveground biomass and LAI of Sasa (an herbaceous species), the inhibition of tree regeneration by Sasa, and the decomposition rates of dead woods and soils. Additionally, the validations of aboveground biomass and tree species composition 32 years after each post-windthrow management effort (dead wood left intact, salvage logging, and scarifications) were conducted. In addition, the dynamics of aboveground biomass and tree species composition in undisturbed forests, the carbon stocks of dead wood and soils, and the net ecosystem production (NEP) were validated in this study for several decades. The calibrations, validations, and the input data of the model are described in detail in Supplementary materials S1, S2, and S3 in Hotta et al. (2021) and Online resource S1 in this study (full of the LANDIS-II input files available in the Supporting data (Hotta 2022)).
Landscape initialization
We focused on the tree species that accounted for 90% of the cumulative biomass in the study landscape and the most dominant herbaceous species as follows: A. sachalinensis, P. jezoensis, P. glehnii, T. japonica, Acer pictum Thunb., B. ermanii, Betula maximowicziana Regel, Quercus crispula Blume var. crispula, Kalopanax septemlobus (Thunb.) Koidz., Fraxinus mandshurica Rupr., Ulmus laciniata (Trautv.) Mayr ex Schwapp., and S. senanensis. Regarding tree species, the initial communities as of 2015, that is, the input data of vegetation, were created based on forest inventory data of UTHF. Regarding herbaceous species, the input data were created by estimating the aboveground biomass of Sasa using the Sasa distribution model developed by Tatsumi and Owari (2013) and tree input data. We input the uniform value 3,100 g m-2 in the study landscape, as this is the mean value of the dead wood data from Hotta et al. (2020), into the initial amount of dead wood. All initial dead woods were assumed to be mostly undecomposed because the data of the amounts of dead woods in each decay class were unavailable. Soil property data were input based on the data from Asahi (1963), who reported soils in UTHF in detail.
Windthrow regime scenarios
Biomass Harvest extension Ver 4.3 (Gustafson et al. 2000) was used to represent various frequencies and intensities of stand-replacing windthrow.
- Historical windthrow regime scenario (Historical)
The scenario in which stand-replacing windthrow occurred with historical frequency (50-year interval) and intensity (20% of the study landscape area) is shown in Table 1. The ratio of windthrow area to the study landscape area of each windthrow event (20%) was determined based on the record of windthrow caused by Typhoon Thad in 1981 (Watanabe et al. 1990). The interval of windthrow (50 years; occurred in 2030 and 2080) was determined based on the data from Abe et al. (2006) (Table 1).
2. Increased windthrow frequency scenario (Frequent)
The scenario in which stand-replacing windthrow occurred with twice the historical frequency (Table 1). The windthrow area ratio to the study landscape area of a single windthrow event was defined as 20% same as the Historical scenario. The interval of windthrow was defined as 25 years (occurring in 2030, 2055, 2080, and 2105).
3. Increased windthrow intensity scenario (Intense)
The scenario in which the stand-replacing windthrow occurred with twice the historical intensity (Table 1). The windthrow area ratio to the study landscape area of a single windthrow event was defined as 40%. The interval of windthrow was defined as 50 years (occurring in 2030 and 2080), which was the same as the Historical scenario.
4. Increased windthrow frequency and intensity scenario (Frequent & Intense)
The scenario in which the stand-replacing windthrow occurred with twice the historical frequency and intensity (Table 1). The windthrow area ratio to the study landscape area of a single windthrow event was defined as 40%. The interval of windthrow was defined as 25 years (occurring in 2030, 2055, 2080, and 2105).
All the live trees within grid cells where windthrow occurred except for advanced seedlings and Sasa were assumed to be destroyed when windthrow occurred because we focused on the stand-replacing windthrow in this study. The destruction ratios of advanced seedlings were dependent on post-windthrow management scenarios, which are described later. Whether the cohorts contained advanced seedlings was determined by the cohort age, and the threshold age was determined by the tree species based on the field data from the windthrow sites in the UTHF caused by Typhoon Thad in 1981 (see Supplementary materials S5 in Hotta et al. (2021)). The grid cells where windthrow occurred were randomly selected within the grid cells where stand age was over 50 years because the windthrow risk increases with the stand age (Everham and Brokaw 1996; Foster 1988) (Online resource S2).
Table 1 The year when windthrow occurs and the ratio of windthrow areas in each windthrow regime scenario
Windthrow regime scenarios
|
2030
(Year 15)
|
2055
(Year 40)
|
2080
(Year 65)
|
2105
(Year 90)
|
Historical
|
20%
|
―
|
20%
|
―
|
Frequent
|
20%
|
20%
|
20%
|
20%
|
Intense
|
40%
|
―
|
40%
|
―
|
Frequent & Intense
|
40%
|
40%
|
40%
|
40%
|
Post-windthrow management scenarios
We considered the following three post-windthrow management practices: (a) Dead woods left undisturbed (windthrow only); (b) salvage logging; and (c) salvage logging and scarification. It is noted that this is the same setting as one used in Hotta et al. (2021)
- Dead woods left undisturbed (windthrow only: WT)
All dead woods generated by the windthrow were left undisturbed, and 20% of advanced seedlings were destroyed due to fallen trees (Table 2). The destruction ratio of advanced seedlings was determined based on the data from the windthrow sites in the UTHF caused by Typhoon Thad in 1981 (Kurahashi et al. 1983).
2. Salvage logging after windthrow (SL)
All dead woods generated by the windthrow were salvaged and removed from forest ecosystems, and 60% of advanced seedlings were destroyed due to salvaging (Table 2). The destruction ratio of advanced seedlings was determined based on the data of the ratio of forest floor area disturbed by logs that were removed (Ohsato et al. 1996).
3. Salvage logging and scarification after windthrow (SLSC)
All dead woods generated by the windthrow were salvaged and removed from the forest ecosystems, and the forest floor in all windthrow areas was scarified after salvaging. All advanced seedlings and 99% of Sasa were destroyed due to scarification (Table 2).
Table 2 Settings of post-windthrow management scenarios (the same setting used in Hotta et al. (2021))
Post-windthrow management scenarios
|
Dead woods generated by the windthrow
|
Advanced seedlings
|
Sasa (dwarf bamboo)
|
Windthrow (WT)
|
Left intact
|
20% destroyed
|
Undestroyed
|
Salvage logging (SL)
|
100% salvaged
|
60% destroyed
|
Undestroyed
|
Salvage logging and scarification (SLSC)
|
100% salvaged
|
100% destroyed
|
99% destroyed
|
Climate scenarios
To consider tree physiological responses to climate change, the current climate and two future climate scenarios were used for simulations. Regarding the current climate scenario (Current), the 1 km mesh climate data (Japanese Meteorological Agency 2012; an average from 1981–2010) in the study area landscape were used for the climate data. We selected the following two extreme scenarios as future climate scenarios: (1) Representative concentration pathway (RCP) 2.6 scenario calculated by the CSIRO-Mk3-6-0 model (hereafter referred to as the RCP2.6 scenario; approximately 1℃ increase in mean annual temperature and almost no change in precipitations by 2100); and (2) the RCP8.5 scenario calculated by the GFDL-CM3 model (hereafter referred to as the RCP8.5 scenario; approximately 6.5℃ increase in mean annual temperature and 10% increase in precipitation by 2100) (Nishimori et al. 2019). In humid climatic zones, temperature increases have a greater impact than precipitation changes. In all of the candidate models considered in this study, precipitation was projected to remain approximately the same or increase in the future. Therefore, it was determined that the impact of drought on forests did not need to be considered in this region. To primarily evaluate the impact of temperature increase, this study adopted the future climate scenarios of RCP2.6 and RCP8.5, which have a small range of variability in precipitation and a large range of variability in temperature. The details of each climate scenario and the reasons for selection are described in Online resource S3.
Simulation and evaluation indicators
We simulated 36 scenarios that were combinations of the four windthrow regime scenarios, three post-windthrow management scenarios, and three climate scenarios. The simulations were replicated three times for each scenario because there were stochasticity and uncertainties related to seed dispersal, cohort establishments, and the selections of windthrow area. The duration of simulations was 115 years from 2015 to 2130.
We evaluated species composition in the forest landscape using aboveground biomass and the number of cohorts successfully established for each species. Net forest-sector carbon balance (Hudiburg et al. 2019; hereafter referred to as NFCB) was calculated to evaluate the difference in the carbon balance between the following two cases: the case when dead woods generated by windthrow were decomposed within forest ecosystems (WT scenario) and the case when dead woods generated by windthrow were salvaged and used as timber and paper products (SL and SLSC scenarios). LCA of management and salvaged woods were conducted to calculate NFCB (Online resource S4). NPP and heterotrophic respiration (Rh) were used to understand the factors that explained the dynamics of the NFCB. Rh is defined as the process of releasing CO2 to the atmosphere by the decomposition of dead wood and soil organic matter. The summary and analysis of the simulation results were conducted in R ver. 3.6.2 (R Core Team 2019).
NFCB is one of the indicators of carbon balance that considers CO2 emissions related to harvesting, transportation, manufacturing, the disposal of the product, and so on (Hudiburg et al. 2019). In this study, we defined the forest sector as shown in Fig. 2. The NFCB was calculated using Equations 1 and 2 after estimating the CO2 emissions related to salvage logging, scarification, the manufacturing of products, and the disposal of products by LCA (Fig. 2, Online resource S4). The CO2 emissions related to the transportation of dead wood were excluded because they were relatively smaller than the CO2 emissions related to the other processes (Owari et al. 2011).
NFCBt = NEPt - CESalvage - CEScarification – CEManufactureLong-lived –
CEManufactureShort-lived - CEManufactureWoodenBoard (t = 0) Eq. 1
NFCBt = NEPt - CEWoodProductsDecomposed t (t ≠ 0) Eq. 2
t: number of years that elapsed since salvage logging; NFCBt: net forest-sector carbon balance in year t; NEPt: net ecosystem production in year t; CESalvage: CO2 emissions related to salvage logging (e.g., use of heavy machinery); CEScarification: CO2 emissions related to scarification (e.g., use of heavy machinery); CEManufactureLong-lived: CO2 emissions related to manufacturing of long-lived products (timber products used for several decades such as building materials and furniture); CEManufactureShort-lived: CO2 emissions related to manufacturing of short-lived products (products to be disposed of in a few years such as wood chips and paper); CEManufactureWoodenBoard: CO2 emissions related to the manufacturing of wooden boards (such as particle boards); and CEWoodProductsDecomposed t: CO2 emissions related to the disposal of products, which were made from dead woods generated by the windthrow, after t years since salvage logging.
The detailed calculations of each CO2 emission are described in Online resource S4. The CO2 emissions related to salvage logging, scarification, and the manufacturing of long- and short-lived products and wooden boards were assumed to occur in the year when windthrow occurred. Regarding the WT scenario, the NFCB was equal to the NEP because all dead woods generated by the windthrow were left intact in forest ecosystems. In this study, scrap woods generated by sawing were used as raw materials for wooden boards to utilize dead woods generated by windthrow as products with the longest possible life. In addition, we analyzed the sensitivity of the NFCB to the usage of scrap wood, and online resource S5 indicated the result in the case when scrap wood was used for paper pulp and chip production, which was the scenario utilizing dead woods generated by windthrow as products with the shortest possible life.