Aboveground biomass increments over 26 years (1993–2019) in an old-growth cool-temperate forest in northern Japan

Assessing long-term changes in the biomass of old-growth forests with consideration of climate effects is essential for understanding forest ecosystem functions under a changing climate. Long-term biomass changes are the result of accumulated short-term changes, which can be affected by endogenous processes such as gap filling in small-scale canopy openings. Here, we used 26 years (1993–2019) of repeated tree census data in an old-growth, cool-temperate, mixed deciduous forest that contains three topographic units (riparian, denuded slope, and terrace) in northern Japan to document decadal changes in aboveground biomass (AGB) and their processes in relation to endogenous processes and climatic factors. AGB increased steadily over the 26 years in all topographic units, but different tree species contributed to the increase among the topographic units. AGB gain within each topographic unit exceeded AGB loss via tree mortality in most of the measurement periods despite substantial temporal variation in AGB loss. At the local scale, variations in AGB gain were partially explained by compensating growth of trees around canopy gaps. Climate affected the local-scale AGB gain: the gain was larger in the measurement periods with higher mean air temperature during the current summer but smaller in those with higher mean air temperature during the previous autumn, synchronously in all topographic units. The influences of decadal summer and autumn warming on AGB growth appeared to be counteracting, suggesting that the observed steady AGB increase in KRRF is not fully explained by the warming. Future studies should consider global and regional environmental factors such as elevated CO2 concentrations and nitrogen deposition, and include cool-temperate forests with a broader temperature range to improve our understanding on biomass accumulation in this type of forests under climate change.


Introduction
Old-growth forests are widely recognized to play an important role in the carbon cycle (Harmon et al. 1990). It has been commonly accepted that old-growth forests are carbon neutral (Odum 1969) and their living biomass is at 'steady state' (Bormann and Likens 1979). However, recent studies indicate that they work as carbon sinks with increasing biomass over centuries (Luyssaert et al. 2008;Tan et al. 2011). Continuous increases in aboveground biomass (AGB) have also been found in temperate (e.g., Foster et al. 2014;Keeton et al. 2011) and tropical (e.g., Baker et al. 2004;Lewis et al. 2009;Phillips et al. 1998;Qie et al. 2017) old-growth forests. Lewis et al. (2004) argued that recent remarkable AGB increases in old-growth Amazonian tropical forests was induced by climate change. On the other hand, in boreal forests biomass growth is more susceptible to climate change in old-growth forests than in young forests, resulting in negative net biomass change (Chen et al. 2016). Thus, assessing climate effects on long-term biomass accumulation in old-growth forests is essential for understanding forest ecosystem functions in a changing climate (McDowell et al. 2020).
Long-term changes in biomass result from the accumulation of short-term changes in the form of gain due to tree growth and loss due to mortality (Hoshizaki et al. 2004). Therefore, to understand how climate affects changes in AGB, the effects of climatic factors on each component need to be taken into account (Chen and Luo 2015;Peña et al. 2018). In addition, endogenous processes such as gap filling in small-scale canopy openings can drive biomass change (McDowell et al. 2020;Phillips et al. 2009): at the local scale, gap formation may cause first a decrease and then an increase in AGB caused by growth promotion of trees around the gap. Repeatedly measured tree census data with tree location can be useful in revealing these processes.
Environmental factors such as topographic position affect both forest biomass (Kubota et al. 2004;Valencia et al. 2009) and tree species composition (Chen and Luo 2015;Kuuluvainen et al. 2017;Ohmann and Spies 1998). For instance, on northern Honshu, Japan, Fagus crenata often dominates forest stands on hillslopes, whereas more tree species occur in riparian areas (Suzuki et al. 2002). Tree species in riparian forests have diverse life history traits (e.g., both shorter and longer lifespans, heavy sprouting (Nakamura and Inahara 2007;Sakio 2020)). Therefore, hillslope and riparian stands are expected to differ in the dynamics (i.e., growth and mortality) and, consequently, the pattern of biomass changes in component species. In addition, a recent analysis of long-term tree census data in northern Japan has revealed different responses among species to changing climate and consequent changes in stand structure and species composition (Hiura et al. 2019). Thus, stands with different topographic characteristics can show different responses to climate change.
Here, we quantify decadal changes in AGB and their processes in relation to endogenous processes and climatic factors, using tree census data measured repeatedly over 26 years  in an old-growth, cool-temperate mixed deciduous forest with different types of topographic units in northern Japan. We ask the following questions: (1) Did AGB show net increase or decrease over the whole forest and study period? (2) Did tree species contribute differently to biomass change among the different types of topographic units? (3) How did gain and loss contribute to the overall changes in stand biomass? (4) Did climatic factors and endogenous processes such as canopy gap formation and recovery influence short-term changes in AGB at the local scale?

Study site
The study was conducted in the Kanumazawa Riparian Research Forest (KRRF) in Yokodake-Maeyama National Forest, Iwate, northern Japan (39° 06′ 37′′ N, 140° 51′ 17′′ E), an old-growth forest with no record of significant anthropogenic disturbances. In KRRF, tree community dynamics have been repeatedly measured since the establishment of a 4.71-ha permanent plot in 1993 ( Fig. 1; Suzuki et al. 2002). This site is one of the core research sites of the Japan Long-Term Ecological Research Network (JaLTER). KRRF has a cool-temperate climate with a mean annual temperature of 8.8 °C and a warmth index (Kira 1991) of 71 °C month (Mahoko Noguchi and Kazuhiko Hoshizaki unpublished data). The mean annual precipitation is 2000 mm, and snow cover lasts 5 months with maximum snow depth of approximately 2 m (Oki et al. 2013). The vegetation depends on the topographic unit. The riparian area is covered with a speciesrich deciduous broadleaved forest consisting of both riparian specialists (Cercidiphyllum japonicum, Aesculus turbinata, Acer mono, Pterocarya rhoifolia, and Ulmus laciniata) and habitat generalists (Fagus crenata and Quercus crispula) (Masaki et al. 2008;Suzuki et al. 2002). The upper slopes and terrace are dominated by F. crenata and Q. crispula. Detailed information on the ecology of component species is available in Hoshizaki et al. (1997Hoshizaki et al. ( , 1999, Masaki et al. (2005), and Osumi (2006). The age of the largest C. japonicum individual is estimated to be more than 500 years (Osumi 2006), indicating that this forest is sufficiently oldgrowth. The natural disturbance regime in KRRF is characterized by canopy gap formation and fluvial sediment movements (Oki et al. 2013). Gap-creating disturbance occurs about every 1-3 years, with gap size ranging from tens to hundreds of square meters (Oki et al. 2013). Recent fluvial sediment movements were recorded in 1988, 1998, and 2007, causing ground disturbance with sizes ranging from 144 to 680 m 2 but no damage to canopy trees.

Field measurement
The 4.71-ha permanent plot was divided to 471 10-m × 10-m quadrats (Fig. 1). The plot ranges in elevation from 400 to 460 m a.s.l., and includes three topographic units: riparian (3.11 ha), terrace (1.03 ha), and denuded slope between them (0.57 ha). In the whole plot, all stems greater than 5 cm in diameter at breast height (DBH) were tagged for identification and mapped, and DBH was measured at the same marked location on each stem in September or October every 2 years from 1993 to 1999 and every 4 years from then to 2019.

Estimation of AGB change and its components
We calculated tree AGB, basal area (BA), and stem density in each topographic unit. Individual tree AGB was estimated by using a general allometric equation for tree species in Japan (Ishihara et al. 2015): (1) ln (y) = −1.196 + 1.162 × ln (D) + 0.338 × ln (D)) 2 − 0.044 × (ln(D) 3 + 0.708 × ln ( ) where y is AGB; D is stem DBH, and ρ is the wood specific gravity of each species (Editorial Board of Wood Industry 1966; Hiroko Kurokawa, Masahiro Aiba and Yusuke Onoda unpublished data; Fujiwara et al. 2007). Confidence intervals of changes in AGB, BA, and stem density were estimated via bootstrapping across 10-m × 10-m quadrats following the method of Valencia et al. (2009). To overview trends in AGB change during the study period and net annual change in AGB, we calculated AGB for three tree size classes: large (≥ 50 cm DBH), medium (15-50 cm DBH), and small (5-15 cm DBH). The net annual change in AGB (in Mg ha −1 y −1 ) was calculated each 4-year period from 1996 to 2019 from the tree DBH data of Fig. 1 Topographic map of the Kanumazawa Riparian Research Forest (KRRF). The solid frame represents the 4.71-ha KRRF plot. Colors denote the three topographic units: blue, riparian (3.11 ha); orange, denuded slope (0.57 ha); green, terrace (1.03 ha). Black dotted lines show the 10-m × 10-m quadrats; thin red lines show the 20-m × 20-m subplots. Contour interval is 2 m 1995-2019. It was then dissected into annual AGB gain and annual AGB loss. Annual AGB gain was calculated separately for growth of trees in each size class and ingrowth, and AGB loss was calculated for mortality of trees in each size class.

Analysis of factors affecting short-term AGB gain at local scale
We examined the effects of climatic condition in each measurement period, canopy gap formation and topography on local-scale AGB gain using two linear mixed-effect models (LMMs). For both models, variables were calculated for every 20-m × 20-m subplot in each 4-year period from 1996 to 2019. Subplot size was determined as an appropriate area to detect gap formation and subsequent recovery in consideration of the range of gap size in KRRF. Each 20-m × 20-m subplot was assigned to one of the three topographic units according to the area of the major topographic unit (based on the number of 10-m × 10-m quadrats therein). When the number of 10-m × 10-m quadrats within a 20-m × 20-m subplot were same between two types of topographic units, we referred to the original map of topographic type and identified the area. Both models used AGB gain as the response variable and included subplot as a random effect.
In model 1, we aimed to investigate whether the amount of AGB gain differed among the measurement periods with the effects of topographic unit and gap formation in the current and previous measurement periods. Fixed effects were initial AGB in the current measurement period, AGB losses in the current and previous measurement periods, topographic unit, and the five 4-year measurement periods between 2000 and 2019, with topographic unit and measurement period as categorical variables. Initial AGB was included as it is expected to be the "capital" for AGB gain by tree growth. AGB losses were indices of gap formation in the current and previous measurement periods.
In model 2, the effect of climate was analyzed separately from the effect of measurement period to avoid multicollinearity. Fixed effects in model 2 were initial AGB in the current measurement period, AGB loss in the current and previous measurement periods, topographic unit, and mean air temperature during the previous autumn (September-November) and the current summer (June-August) over the measurement periods. For example, mean air temperature during the previous autumn of the measurement period 1996-1999 is the mean air temperature during September-November, 1995, and that of the current summer is the mean during June-August, 1996-1999.Both types of mean air temperature have a major influence on annual DBH growth of individual trees in most dominant species of KRRF (Michinari Matsushita, Daiki Sugiura and Kazuhiko Hoshizaki manuscript in preparation). As the on-site temperature data do not cover the entire study period, we used data from the nearest weather station, at Wakayanagi (39° 08′ N, 141° 04′ E; 97 m a.s.l.: Japan Meteorological Agency, https:// www. data. jma. go. jp/ gmd/ risk/ obsdl/ index. php), 18 km east of the study site. In addition to the aforementioned factors, interactions between the topographic unit and climatic factors were included in the model's fixed effects of the model to examine whether different topographic characteristics shows different responses to climatic factors (model 2.1).
These models were fitted by the lme4 v.1.1-21 package (Bates et al. 2015) in R 3.6.3 (R Core Team 2020). All variables except for categorical variables (i.e., topographic unit and measurement period) were standardized before the analyses. To evaluate the variance explained by the models, we calculated two R 2 values for mixed-effect models following the method of Nakagawa and Schielzeth (2013) and Nakagawa et al. (2017): marginal R 2 (R 2 LMM(m) ), which is the proportion of the variance explained by fixed effects, and conditional R 2 (R 2 LMM(c) ), which is the proportion of the variance explained by both fixed and random effects. These were calculated by the MuMIn v. 1.43.15 package (Bartoń 2020) in R.

Overall changes in AGB at plot scale
BA in 1993 was greatest in the riparian unit (34.2 m 2 ha −1 ) and least in the denuded slope unit (Table 1). From 1993 to 2019, BA increased significantly in all topographic units. AGB was greatest in the terrace unit (246.0 Mg ha −1 ) at the beginning of the study period (Table 1). It increased significantly in all topographic units during the study period, increasing in most 4-year periods except for some short pauses; for instance, from 2011 to 2015 in the riparian and denuded slope units (Fig. 2). AGB of large trees (≥ 50 cm DBH) in 1993 occupied 76.7% of total AGB in riparian, 70.6% in denuded slope, and 77.7% in terrace units. Trends of increasing total AGB in the riparian and terrace units were similar to those of large-tree AGB. During the study period, stem density declined in the riparian and terrace units but increased in the denuded slope unit (Table 1). The change in stem density was significant only in the riparian unit.
In the riparian unit, C. japonicum had the largest AGB at the beginning of the study period, followed by F. crenata, A. turbinata, Q. crispula, and A. mono ( Table 2). AGB of these species, except for Q. crispula, increased during the study period. Pterocarya rhoifolia had the greatest increment in AGB over the study period, accounting for 52.3% of the total increment in the riparian unit, followed by A. turbinata at 25.4%. In contrast, several other species with relatively small AGB at the beginning, such as Zelkova serrata and Ulmus laciniata, showed a decline in AGB during the study period. The denuded slope and terrace units were dominated by F. crenata and Q. crispula, and the denuded slope by A. mono as well (Table 2). All these species had an increase in AGB during the study period, maintaining the AGB-based rank of species composition.
Annual gain in AGB remained at approximately 3 Mg ha −1 year −1 with some differences among the measurement periods: larger in 2008-2011 and 2016-2019 in all topographic units (Fig. 3). In the riparian and terrace units, large-and medium-sized trees accounted for most of the annual gain. Annual losses in AGB fluctuated among the 4-year periods, and were largest in 2012-2015 in all topographic units. Regardless of topographic unit, measurement periods with greater loss of AGB of large trees tended to have greater total loss of AGB. As a consequence, net annual change in AGB ranged from − 0.6 to + 2.6 Mg ha −1 year −1 but stayed positive except in 2012-2015 in the riparian and denuded slope units.

Effects of climate and gap formation on short-term AGB gain at local scale
Local-scale AGB gain in the 20-m × 20-m subplots was positively influenced by initial AGB in each measurement period with the largest effect size among the numeric explanatory variables (LMM model 1: Table 3, Fig. 4). It was significantly greater in subplots with larger AGB loss in the previous measurement period but smaller in subplots with larger AGB loss in the current measurement period. It was not significantly affected by topographic unit. It differed significantly among the measurement periods: smaller in 2004-2007 and larger in 2008-2011 and 2016-2019. In model 1, R 2 LMM(m) = 0.31 and R 2 LMM(c) = 0.76, indicating that 31% of the variation was explained by fixed effects and 76% by fixed and random effects. In model 2 (Table 4), the effects of initial AGB, AGB loss in the current and previous measurement periods, and topographic unit were almost identical to those in model 1. Local-scale AGB gain was larger in measurement periods with higher mean air temperature during the current summer but smaller in those with higher mean air temperature during the previous autumn. The absolute effect size of these two variables was almost equivalent. Model 2 explained almost identical variation as model 1, with R 2 LMM(m) = 0.30 and R 2 LMM(c) = 0.75. None of the interactions between the topographic unit and climatic factors (i.e., mean air temperature during the previous autumn and the current summer) were significant (model 2.1, Table S2). Additionally, the interactions did not improve performance of the model with R 2 LMM(m) = 0.30 and R 2 LMM(c) = 0.75 in model 2.1.

Discussion
AGB of KRRF increased steadily over the 26 years in all topographic units, with increments of 30-35 Mg ha −1 in each (Table 1). BA also increased over the study period, even though it was initially equivalent to values reported in other cool-temperate old-growth forests in Japan (Masaki et al. 1992;Nakashizuka 1988;Seiwa et al. 2013), indicating that the forest had already been well stocked. These results are consistent with reports that temperate old-growth forests continuously gain biomass over the long term (Keeton et al. 2011;Luyssaert et al. 2008). This continuous stand-scale biomass increment was attributable mainly to an increase in AGB of large trees, in agreement with the reported global importance of large trees in determining stand AGB (Lutz et al. 2018;Slik et al. 2013). Patterns of tree growth or stand-biomass-change vary across tree species composition and diversity, as well as with environmental conditions such as topography (Kubota et al. 2004;Valencia et al. 2009). In KRRF, topography has been reported to determine tree species distribution through affecting seedling survival and growth differently across species (Masaki et al. 2005). They argued that decreased seedling survival and growth in the terrace unit was caused by lower water availability. However, a steady increase in AGB was common to all three topographic units (Table 1; Fig. 2). Furthermore, neither topographic units (Tables 3 and  4) nor their interaction with climate (Table S2)   influence on local-scale AGB gain. These results suggest that topography-induced differences in water availability do not contribute to spatio-temporal variations in AGB gain of adult trees in KRRF under the recent climate. We attribute the synchronous increase in AGB among the topographic units with different tree species composition primarily to the AGB increment in F. crenata, which is dominant in all three topographic units (Table 2). A growing abundance of F. crenata has been documented in several stable old-growth forests (Seiwa et al. 2013;Yamamoto and Nishimura 1999). Increases in both AGB and stem density of F. crenata in KRRF may be due to lack of remarkable disturbance even in the riparian unit during the study period. In the riparian unit, however, the contribution of F. crenata to the AGB increment is lower than in other topographic units, partially due to the smaller initial dominance of this species. Instead, Pterocarya rhoifolia, a riparian specialist of cool-temperate forests in Japan (Sakio et al. 2002), made the largest contribution to the stand AGB increment (Table 2). Despite the substantial decline in its stem density (Table S1), its AGB at the end of the study period was 3 times the initial value. It is likely that the fast growth of P. rhoifolia (Sakio 1993) is associated with its rapid increase in AGB. Although AGB decreased in some species such as Z. serrata and U. laciniata in the riparian unit, P. rhoifolia compensated for the decrease and resulted in the stand-level AGB increase.
Mortality is the major cause of reduced growth or decline in AGB (Schuster et al. 2008;Xu et al. 2012). Although large disturbances such as strong typhoons, insect outbreaks, or severe flooding in the riparian unit were not recorded during the 26 years, the loss of AGB in KRRF varied substantially among the topographic units and among the 4-year measurement periods (Fig. 3). These variations were explained mainly by the spatio-temporal variation in mortality of large trees. A significant contribution of large-tree mortality to the AGB loss has also been reported in other old-growth forests (Hoshizaki et al. 2004). A large amount of AGB loss in 2012-2015 is due to the mortality of largersized trees in this period (Mahoko Noguchi and Kazuhiko Hoshizaki unpublished data). Despite these temporal and spatial variations, the AGB loss generally remained smaller than the AGB gain, bringing about a positive change in AGB in most of the measurement periods. The temporal change of stand-level AGB appears to be inconsistent with the assumed  long-term balance between biomass loss caused by canopy gap formation and subsequent gain during gap recovery. In contrast to AGB loss, temporal fluctuations in AGB gain tended to synchronize across the topographic units at the stand scale (Fig. 3). The results of model 1 indicate that local-scale AGB gain also differed among measurement periods even after adjustment for initial AGB and disturbance during each period (Table 3). As expected, initial AGB positively influenced local-scale AGB gain (Fig. 4). Larger AGB loss in the previous measurement period caused greater AGB gain, suggesting that variations in local-scale AGB gain are partially explained by recovery in and around canopy gaps. Local-scale AGB gain also substantially differed among the measurement periods. The results of model 2 suggest that the observed temporal variations in AGB gain are caused by climatic factors: a warmer current summer had positive effects whereas a warmer previous autumn had negative effects on AGB gain (Table 4). The positive response of tree growth rates to high temperature in the growing season has been reported in deciduous broadleaved species in cool-temperate forests in northern Japan (Hiura et al. 2019). In KRRF, this response is consistent with that of individualbased annual tree growth (Matsushita et al. manuscript in preparation), which often shows considerable inter-annual variations (Ohtsuka et al. 2009) in response to climatic factors (Nabeshima et al. 2010). The negative effect of warmer previous autumn on annual tree growth of deciduous broadleaved species is also found in KRRF (Matsushita et al. manuscript in preparation) and in a cool-temperate forest in central Japan (Shen et al. 2020). Also, over the Northern Hemisphere, autumn warming causes a larger increment in respiration than in photosynthesis and leads to net carbon loss (Piao et al. 2008).
Our models incorporating mean air temperature explained a considerable amount of variation in local AGB gain, although the analysis did not include other potential factors that enhance tree growth such as change in precipitation (Hiura et al. 2019). Both the summer and autumn temperatures at the weather station nearest to KRRF have shown a substantial rise over the past 40 years (Fig. S1). As the positive and negative effects of summer and autumn temperatures were equivalent in absolute size (Table 4), the influences of summer and autumn warming on AGB growth appeared to be counteracting. Therefore, the observed steady AGB increase in KRRF is not fully explained by decadal trends of warming. The simulation results of European temperate forests show that elevated CO 2 concentrations and nitrogen deposition underpins current high stand growth under the recent climate (Pretzsch et al. 2014). Future studies should consider these global and regional environmental factors, and include cool-temperate forests with a broader temperature range to improve our understanding on biomass accumulation in this type of forests under climate change.