The impact of climate variation and disturbances on the carbon balance of forests in Hokkaido , Japan

Introduction Conclusions References


Introduction
Changes in carbon flux and storage in forest ecosystems are influenced by climate at various temporal and spatial scales.However, they are also affected heterogeneously by artificial and natural disturbances at the local scale.The time course of net ecosystem production (NEP) is affected by climatic conditions, because biochemical processes such as photosynthesis, respiration, phenology, and allocation all respond to climate (Hirata et al., 2008;Saigusa et al., 2008;Reichstein et al., 2013).Disturbance events also drastically change NEP and carbon storage.Forests in boreal and tropical regions frequently lose large carbon stocks through wildfires (Pan et al., 2011;Hirano et al., 2012;Ueyama et al., 2013).Damage by insects has had a similar impact on forest ecosystems (Kurz et al., 2008a).Several years after disturbance, however, regenerated forests become carbon sinks because the NEP of a young forest is higher than that of an old forest (Pan et al., 2011).After several decades, NEP tends to approach an equilibrium state (Janisch and Harmon, 2002;Goulden et al., 2011).The drastic changes of carbon flux caused by disturbance in forest ecosystems affect global climate (Dale et al., 2001;Turner, 2010).Although large carbon emissions caused by disturbance may accelerate global warming, enhanced carbon sinks after disturbance mitigate global warming (Pan et al., 2011).Despite these previous findings, the effects of forest disturbance on global climate remain to be clarified (Kurz et al., 2008b).
The effects of disturbance and subsequent recovery on the carbon balance are attributable to carbon uptake by regrowing trees and carbon loss by decomposition of dead wood (i.e., residues; Janisch and Harmon, 2002;Grant et al., 2010;Goulden et al., 2011).Therefore, controlling the degree of disturbance and recovery of forest vegetation and managing residues are expected to mitigate global warming by reducing carbon emission and enhancing carbon sequestration (Liu and Han, 2009).By quantifying the mechanisms underlying the carbon balance in forest ecosystems and the effects of carbon management practices, researchers can apply these lessons to the mitigation of climate change (as outlined in the Kyoto Protocol, Clean Development Mechanism, and Reducing Emissions from Deforestation and Forest Degradation), improve conservation, devise sustainable management techniques for forests, and enhance forest carbon stocks in developing countries (REDD+).From a long-term perspective, it is easier to maintain or enhance the function of carbon sinks in managed forests than in unmanaged forests, thus the management of carbon sinks was adopted in the Kyoto Protocol and the subsequent international framework under the Durban Platform.
The eddy covariance technique is a useful tool for investigating the interaction between terrestrial ecosystems and the atmosphere and evaluating the carbon balance.Despite the technique's utility, understanding these systems requires the continuous and long-term measurement of fluxes.The longest-running measurement of eddy flux was started in 1990 in Harvard Forest (Urbanski et al., 2007), and the measurement periods at most sites in Asia are less than 10 years (Mizoguchi et al., 2008).In addition, eddy flux measurements just after disturbance have been conducted at only a few forest sites (Kowalski et al., 2004;Takagi et al., 2009).However, understanding the effect of disturbance on the carbon balance requires data gathered over several decades, from the time the disturbance occurs until the forests reach maturity (Janisch and Harmon, 2002;Grant et al., 2010;Goulden et al., 2011).The changes in carbon stock and carbon flux caused by disturbance have been investigated by using the chronosequence approach, in which flux measurements are made in similar forest stands that differ in age (Howard et al., 2004;Amiro et al., 2006Amiro et al., , 2010;;Zha et al., 2009;Goulden et al., 2011).Numerical simulation using ecosystem models also has been used to evaluate the longterm effects of disturbance and recovery on the carbon balance (Janisch and Harmon, 2002;Taylor et al., 2008;Ito, 2012).Studies combining the chronosequence approach and ecosystem models have been conducted (Law et al., 2004;Grant et al., 2010).
Our objective was to evaluate the long-term (about 50year) effect of climate, disturbance, and subsequent recovery on the carbon balance of forest ecosystems.We hypothesize that (1) the carbon balance of a forest ecosystem is controlled by the amount of residues, (i.e., the dead biomass produced by the disturbance and left at the site); and (2) the period required for the carbon balance to be restored can be controlled by managing these residues.The study sites were artificial larch (Larix spp.) forests in Hokkaido, Japan, that had been planted after the clear-cutting of mixed forest stands.First, we validated the ecosystem model by using flux data obtained from stands of a wide range of ages, including an old forest before disturbance, a young forest just after disturbance, and a middle-aged forest.Second, we used the model to investigate which factors influenced the long-term carbon balance from the young forest stage to the middle-aged stage.Finally, we analyzed the sensitivity of the carbon balance to the amount of residues.

Site description
Carbon flux data were obtained from two cool temperate larch forests at different ages and a mature mixed forest in Hokkaido, which is the northern main island of Japan.One was a young larch forest planted after the clear-cutting of a conifer-hardwood mixed forest in the Teshio Experimental Forest of Hokkaido University (TSE site; 45 • 03 N, 142 • 06 E), and the other was a middle-aged larch forest at the Tomakomai flux research site (TMK site; 42 • 44 N, 141 • 31 E).The characteristics of the Teshio and Tomakomai sites are listed in Table 1.
The Gleyic Cambisol soil at the Teshio site is more fertile than the volcanogenic Regosol soil at the Tomakomai site.At both sites, the forest floor is covered with a thick understory, which consists mainly of evergreen dwarf bamboo at Teshio and ferns at Tomakomai.Details of the Teshio and Tomakomai sites were given by Takagi et al. (2009) and Hirata et al. (2007), respectively.
The Teshio site had been a mature (∼ 200-year-old) mixed forest in the late successional stage.In December 1972, it was damaged by strong wind in a blizzard and some trees were blown down.From January to March 2003, the trees at the site were clear-cut, with 55 % of the logs removed and the residual biomass remaining in the forest.The Sasa understory was then strip-cut into 4 m wide strips, and 2-yearold hybrid larch trees were planted in the spaces between the strips (Takagi et al., 2009).At the Tomakomai site, Japanese larch was planted in 1958 after the forest had been devastated

Data processing
In the model simulations (for model description, see Appendix A), we used CO 2 flux data collected by using the eddy covariance method at flux towers, as well as climate data gathered at the Teshio and Tomakomai sites.Details of the tower equipment and instruments are given by Takagi et al. (2009) for Teshio and by Hirata et al. (2007) for Tomakomai.The eddy covariance method provided half-hourly NEP measurements, which were partitioned into gross primary production (GPP) and ecosystem respiration (RE).RE was obtained from the relationship between nighttime NEP and temperature, and GPP was obtained by subtracting RE from NEP.The data processing, quality control, gap-filling, and partitioning into GPP and RE are described by Hirata et al. (2008) and Aguilos et al. (2014).We used 1 year of data (2002) for the ∼ 200-year-old mixed forest at the Teshio site; 10 years of data (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) for the young larch forest at Teshio, which ranged in age from 1 to 10 years old during that period; and 3 years of data (2001)(2002)(2003) for the middleaged larch forest at the Tomakomai site, which ranged in age from 44 to 46 years (Table 1).

Model application
To simulate the carbon cycles of the larch forests at Teshio and Tomakomai, we adapted the Vegetation Integrative Simulator for Trace Gases (VISIT) model (Ito, 2010), which was developed on the basis of a simple carbon cycle model (Sim-CYCLE; Ito and Oikawa, 2002) and subsequently improved (Ito et al., 2005).Details of VISIT are described in Appendix A.
The VISIT model has been modified by using flux tower data obtained from the AsiaFlux database, including the middle-aged larch forest at the Tomakomai site and mature mixed forest at the Teshio site, and carbon fluxes simulated with the amended model have shown good agreement with the observed values (Ito, 2008;Ichii et al., 2010Ichii et al., , 2013)).However, the VISIT model has not been validated for young forests during the initial stage of growth.Takagi et al. (2009) reported that the understory played an important role in the carbon balance of the young larch forest at the TSE site.Under the scheme of the conventional VISIT model (Ito et al., 2006), however, the understory was treated as evergreen vegetation and did not show apparent seasonal variations.Before we conducted the analysis, we modified the model to amend the shortcomings noted by Takagi et al. (2009).To strengthen the role of the understory in the carbon balance, we modified the phenology of the understory by treating it as deciduous plants and changing the parameter values for the specific litter fall rate (Eq.A6).As a result, in the modified VISIT model, the leaf area index (LAI) of the understory R. Hirata et al.: Impact of climate variation and disturbances showed seasonal variation and the understory absorbed more CO 2 and grew faster than in the conventional model.
To simulate long-term carbon fluxes at the Tomakomai site, we used long-term daily meteorological data derived from reanalysis of a global climate data set produced from 1948 to 2010 after merging of the US National Centers for Environmental Prediction and the US National Center for Atmospheric Research (NCEP/NCAR; Kistler et al., 2001) and Climatic Research Unit time-series data sets (CRU TS3.1; Mitchell and Jones, 2005) corrected by in situ observational data.Details of the climate data set are given by Ichii et al. (2013).The model was initialized by a spin-up run for 2000 years to create the dynamic equilibrium of soil organic matter and vegetation components using the reanalysis climate data corrected by in situ observational data from 1948 to 2010 for the Teshio and Tomakomai sites.CO 2 concentrations were constant at 312 ppm until 1948.After the spin-up run, a simulation was conducted for the observation period (11 years for Teshio and 3 years for Tomakomai sites) and for the long-term period for Tomakomai.The annual atmospheric CO 2 concentration was linearly increased from 312 ppm in 1948 to 404 ppm in 2012 for Teshio, from 312 ppm in 1948 to 380 ppm in 2003 for Tomakomai, and from 312 ppm in 1948 to 398 ppm in 2010 for the 52-year simulation for Tomakomai.We implemented disturbance events that matched those in the sites' histories as much as possible.In the case of the Teshio site, the mixed forest was devastated in 1800 and regrew: 90 % of the foliage and stems and 80 % of the roots were exported, and residues were moved to the litter compartments as dead vegetation.Because of the blizzard in 1972, half the trees were removed.To simulate the conversion of mixed forest to larch forest (clear-cutting and planting larch) in 2003, 55 % of the stems were removed as wood harvest and the residues remained on the forest floor as dead stems (Takagi et al., 2009).Whole foliage and roots were left on and in the forest floor, respectively (Takagi et al., 2009).In the model, live biomass of foliage, stems, and roots is reduced to zero by clear-cutting, and residues are treated as carbon in litter compartments that are moved from the live biomass.Therefore, 55 % of the live stem biomass was removed and treated as emission, and 45 % of the live biomass was moved to the stem litter compartments.All live foliage and roots were moved to the foliage litter compartment and the root litter compartment, respectively.When a conversion event occurred, the parameters of the mixed forest simultaneously changed to those of larch forest.In the case of Tomakomai, the information before the larch plantation was initiated in 1958 and the details of thinning were not clear.Therefore, we applied the same conditions as at the Teshio site to the Tomakomai site before plantation initiation.The mixed forest was devastated in 1760 and regrown.We implemented a disturbance in 1930, and half the trees were exported.Mixed forest was converted to larch forest (clear-cutting and planting larch) in 1958.When thinning was implemented in 1970, 1985, and 2004, 30 % of the live stems were exported and 30 % of the live roots were moved to the litter compartments in the model.In the case of larch forest, which is a deciduous needle-leaf forest, no live foliage was moved to the litter compartments because thinning was conducted during dormancy.The intensity of thinning followed the common practice in Japan.

Management and climate scenarios
A series of scenarios was designed to examine the attribution of long-term carbon fluxes to disturbance and individual climate factors at the Tomakomai site.We compared carbon fluxes calculated by using historical climate data with constant scenarios that had no interannual variation in each climate factor.For each scenario, a series of individual factors was replaced by an ensemble average for 1948-2010 and the same seasonal variation was retained during the simulation period.Constant scenarios started in 1958, when clearcutting occurred and the plantation was initiated.We focused on three types of disturbance: (1) clear-cutting; (2) conversion, in which the mixed forest was clear-cut and larch was planted; and (3) thinning, in which 30 % of the forest biomass was cut and removed.VISIT was run to examine the effect of disturbance under 10 scenarios: full scenario (S full ), in which temperature, precipitation, solar radiation, vapor pressure deficit (VPD), and CO 2 concentration are historical with historical disturbances such as clear-cutting, conversion and thinning, climate constant (S const−climate ), in which temperature, precipitation, solar radiation, and VPD are constant; temperature constant (S const−Ta ); precipitation constant (S const−precipitation ); solar radiation constant (S const−Sd ); VPD constant (S const−VPD ); and CO 2 constant (S const−CO2 ).These scenarios were run with clear-cutting and conversion of mixed forest to larch forest and without thinning.We also conducted a non-conversion scenario (S non−conv ), in which a mixed forest without clear-cutting continued to exist (neither clear-cutting nor plantation occurred); a nonclear-cutting scenario (S non−cut ), in which a larch forest without clear-cutting continued to exist; and a non-thinning scenario (S non−thin ) (Table S1 in the Supplement).The effects of climate, temperature, precipitation, solar radiation, VPD, CO 2 , conversion, and clear-cutting were evaluated by subtracting S const−climate , S const−Ta , S const−precipitation , S const−Sd , S const−VPD , S const−CO2 , S non−conv , and S non−cut , respectively, from S non−thin .The effect of thinning was evaluated by subtracting S non−thin from S full , respectively.

Ecosystem compensation points
To evaluate the recovery of carbon fluxes after clear-cutting disturbance, we used three indices of ecosystem compensation points.These indices provide useful information for carbon management related to disturbance.
To understand when the ecosystem turns into a carbon sink in the recovery phase on an annual basis, we used the ecosystem compensation point for annual NEP (ECP NEP ), which is the number of years until the ecosystem shifts from being an annual carbon source as a result of disturbance to being an annual carbon sink.Kowalski et al. (2004) first used this term, but Thornton et al. (2002) called it the "length of the source period." Despite the ecosystem reverting to net CO 2 fixation at the ECP NEP , the ecosystem was still under a carbon debt because carbon was released after the disturbance before the ECP NEP .Therefore, we need to account for the carbon release before the ECP NEP when we consider carbon management after disturbance.Here, we use the ecosystem compensation point for cumulative NEP (ECP CNEP ), which is the number of years until cumulative NEP shifts from negative to positive after the disturbance; this represents how long it takes to pay back the carbon released after the disturbance.
If cut logs or fallen trees burn and the carbon is released into the atmosphere immediately, we also need to take this carbon into account.The third index is the ecosystem compensation point for cumulative net ecosystem carbon balance (NECB) (ECP CNECB ), which is the number of years until cumulative NECB shifts from negative to positive after the disturbance.NECB is calculated by subtracting removed carbon from NEP and adding imported carbon.

Sensitivity test
We conducted a sensitivity test by using long-term data (from 1958 to 2010) of a larch forest in the Tomakomai site for the strength of disturbance.The proportions of the carbon emission were set to three scenarios: (1) emission of tree biomass is 0 % and that of residue is 100 % (S E0R100 ); (2) emission of tree biomass is 50 % and that of residues is 50 % (S E50R50 ); and (3) emission of tree biomass is 100 % and that of residues is 0 % (S E100R0 ).To simplify the test and the results, emission of tree biomass included foliage, stems, and roots, although this was unrealistic because only the stem is exported and most of the foliage and stems remain on the forest floor in the actual forest management, as mentioned in Sect.2.3.Residues of foliage, stems, and roots were moved to each litter compartment.

Model validation
The results of the simulation were compared with NEP observed by the eddy covariance technique and with GPP and RE, which were estimated from the observed NEP by using semi-empirical models (Fig. 1, Table S2 in the Supplement).The mature mixed forest absorbed CO 2 during the growing season, but NEP shifted from positive to negative after clearcutting and planting of larch trees in 2003 because of an abrupt decrease in GPP (Fig. 1).After that, CO 2 release from the ecosystem decreased with the increase in GPP.Field ob- servations showed that apparent CO 2 absorption during the growing season began in 2008, and the ecosystem reverted to fixing CO 2 annually (0.5 t C ha −1 year −1 ) in 2010.On the other hand, the simulated annual NEP changed from negative to positive in 2009 (1.6 t C ha −1 year −1 ).GPP showed a drastic change after clear-cutting, whereas the change in RE was small according to both observations and the simulation (TSE in Fig. 1).At the Tomakomai site, the peak of NEP occurred in early summer and then a rapid decrease was observed in midsummer (Hirata et al., 2007), and the model values also followed this pattern (TMK in Fig. 1).We performed statistical tests for comparisons of carbon fluxes calculated from all three forest types-the mature mixed forest (2002 in Teshio), young larch forest (2003( -2012 in Teshio) in Teshio), and middle-aged larch forest (2001)(2002)(2003) in Tomakomai)-with those estimated from field observations on daily, monthly, and yearly scales (Table S2 in the Supplement).However, in the case of yearly statistics, only NEP and GPP are shown for the young larch forest, because too few data were available or there were no significant relationships for other fluxes.Although daily statistics showed some biases (NEP of the three forest types, RE of mature mixed forest, for which the slopes were less than 0.8; Table S2 in the Supplement) or weak relationships between simulated and observed carbon fluxes (NEP of the mature mixed forest and the young larch forest, for which the R 2 values were less than 0.5; Table S2 in the Supplement), most of the monthly model output of carbon fluxes agreed well with observed values (the slopes were over 0.8 or less than 1.2 and the R 2 values were over 0.5: Table S2 in the Supplement).In this model, using monthly values appeared to be most suitable for evaluating seasonal variations in carbon fluxes, rather than using daily values.On the whole, the high correlation of the monthly data indicated that the model successfully simulated the seasonal patterns of observed NEP, GPP, and RE.
On the yearly scale, model outputs of annual NEP and GPP of the young larch forest matched the observed values well (R 2 values were 0.82 and 0.62 for NEP and GPP, respectively; Table S2 in the Supplement), although there was some bias in NEP (slope was 1.51 for NEP; Table S2 in the Supplement).In contrast, we did not obtain significant correlations of annual NEP, GPP, or RE for mature mixed forest and middle-aged larch forest or for RE of young larch forest between the model estimates and the observations, because the differences among annual values were small (Table 2).However, the average of each annual carbon flux estimated by the model was similar to the observed value (Table 2).Overall, the model output was consistent with the observed fluxes from the seasonal to the yearly scale.We therefore used the model to evaluate the long-term effects of disturbance on carbon fluxes.
The simulated biomass of trees (60 t C ha −1 ) in the middle-aged larch forest at Tomakomai was slightly larger than the observed value (52 t C ha −1 ).The simulated plant biomass of the mature mixed forest (94 t C ha −1 ) was smaller than the observed value (99 t C ha −1 ; Aguilos et al., 2014).The simulated vegetation carbon increment of the young larch forest at Teshio from 2003 to 2009 was 22 t C ha −1 ; this was slightly larger than the observed value (17 t C ha −1 ; Aguilos et al., 2014).

Management and climate scenarios
We examined the simulated long-term NEP, from 1948 to 2010 (Fig. 2a) at the Tomakomai site.There was an abrupt shift from carbon sink to source after clear-cutting and conversion of the mixed forest to a larch plantation, followed by a decrease in carbon release, and then a return to being a carbon sink.This pattern was similar to the interannual pattern of NEP observed in the young larch forest at Teshio from 2002 to 2012.After that, NEP increased asymptotically, with fluctuations.

Effect
First (1958-1967) Second (1968-1977) Third (1978-1987) Fourth (1988-1997)  We analyzed how carbon fluxes were affected by climate, CO 2 , clear-cutting, conversion, and thinning (Fig. 2, Table 3).The effect of disturbance and climate factors on annual NEP, as well as the long-term trends (i.e., over the 52-year period from 1958 to 2010) of each effect were evaluated by time-averaging the effect for each decade after clear-cutting in 1958 (Table 3).However, by calculating the time-average, the positive and negative fluctuations of the effects compensated for each other; therefore, the standard deviation of the effect was also used to evaluate the effect of interannual variation in each factor (Table 3).Annual NEP declined immediately after clear-cutting (Fig. 2a), and the effect of this disturbance was much larger than the effect of climate variation (Fig. 2b, Table 3).In the second decade after clear-cutting, the effects of clear-cutting on annual NEP shifted from negative to positive (Fig. 2b, Table 3).Clear-cutting still affected the average NEP, and the effect was larger than that of climate even 50 years after this disturbance (Fig. 2b, Table 3).However, climate had a relatively large effect on the interannual variation of NEP 30 years after clear-cutting (Table 3).This indicates that the year-to-year variation in the annual NEP was driven by climate variation, whereas the 52-year trend in NEP resulted from the disturbance that had occurred 52 years previously.
The effect of clear-cutting on annual GPP approached zero about 15 years after the disturbance (Fig. 2c), suggesting that after this period the GPP of the post-harvest forest was similar to that of intact forest.The effect of clear-cutting on annual RE was continuously negative, even after 50 years (Fig. 2d).Even though the GPP (or net primary production, NPP) of clear-cut forest was similar to that of intact forest, there was less dead biomass, which determines heterotrophic respiration (RH; equation A5), in the disturbed forest.Consequently, clear-cutting reduced RH (Fig. 2e).Thus, although GPP (or NPP) recovered to a level similar to that of intact forest 52 years after clear-cutting, the decomposition rate was reduced because the amount of litter decreased.
The effect of conversion to a plantation on all annual carbon fluxes was positive a decade after clear-cutting (Fig. 2be, Table 3).Conversion to a larch plantation had a large effect on annual GPP, RE, and RH (Fig. 2c-e), because GPP and RH in the larch forest were larger than the corresponding values in a mixed forest (S non−thin vs. S non−cut after disturbance: 1031 vs. 873 t C ha −1 for GPP, 529 vs. 428 t C ha −1 for RH).However, the effect of conversion on annual NEP was less than those on other flux components (Fig. 2b) because increased GPP was offset by increased RE.
Although the first thinning (1970) decreased annual NEP, the second (1985) and third thinning ( 2004) enhanced annual NEP 3 years and 1 year after thinning, respectively (Fig. 2b, Table 3).The reason for enhanced NEP is that thinning events decreased RE and RH in 1985 and 2004 (Fig. 2d,  e), whereas the effect on GPP was small (Fig. 2c).The reason for decreased RH is that soil carbon of S full was smaller than that of S non−thin .In the model, when thinning occurred, soil carbon increased by moving root biomass to the litter compartment.However, root litter is easily decomposable (decomposition rate 1 year after thinning were 36, 25, and 28 % in 1970, 1985 and 2004, respectively).In addition, litter supply from live trees decreased by thinning.Therefore, soil carbon of S full was smaller than that of S non−thin 3 years (1985) or 1 year (2004) after thinning.The reason why GPP did not change drastically is that the extra growth of the remaining trees compensated for the removed trees.
The anomaly of the effect of climate can be separated into each meteorological component (Table 3).The abrupt decrease in NEP in 1984 (Fig. 2b) was caused by the fact that there was much less precipitation than normal (694 mm vs. the average of 1000 mm).At this time, GPP also abruptly decreased, whereas RH did not change.Fifty years after clearcutting, the standard deviation of NEP, which is an indicator of interannual variation in NEP, was affected mainly by air temperature (Table 3).The effects of CO 2 on annual productivity (NEP and GPP) and respiration (RE and RH)  and (c) cumulative net ecosystem carbon balance (NECB) at the Tomakomai site.These were simulated on the basis of the four treatment scenarios for clear-cutting and residues.Clear-cutting was not implemented in the non-disturbance scenario.S E0R100 : emission of logs 0 %, that of residues 100 %; S E50R50 : emission of logs 50 %, that of residues 50 %; S E100R0 : emission of logs 100 %, that of residues 0 %.

Ecosystem compensation points
The ecosystem shifted from a net carbon source to a sink 5 years (ECP NEP ) after clear-cutting at Tomakomai site (Fig. 2a).It took 14 years for cumulative NEP to shift from negative to positive after the disturbance (i.e., ECP CNEP ).It took 23 years for cumulative NECB (subtracting emission as logs from NEB) to shift from negative to positive after the disturbance (ECP CNECB ).

Sensitivity test
We examined the sensitivity of carbon flux to the residual biomass at the Tomakomai site after clear-cutting. Figure 3 shows the results of the sensitivity analysis of carbon balance in terms of changes in the proportions of residual debris as foliage, stems, and roots.Short-term effects are given in  The second row shows the percentages relative to that in the control experiment (S E100R0 ) as 100 %.The third row shows the cumulative values of NEP, GPP, RE, and RH from the clear-cutting in 1958 to 2010.The fourth row shows the percentages relative to that in the control experiment (S E100R0 ) as 100 %.When all of the logs were left on the forest floor (S E0R100 ), the ecosystem released more carbon immediately after clearcutting than in the other two scenarios (S E50R50 and S E100R0 ), and the least carbon was released when all the logs were exported (S E100R0 ; Fig. 3a, NEP of S E100R0 in Table 4).Therefore, the longest and the shortest ECP NEP were found in scenarios S E0R100 (8 years) and S E100R0 (4 years), respectively (Fig. 3a, ECP NEP in Table 4).ECP NEP of S E50R50 had an intermediate value of 5 years.
Over the 52 years, cumulative NEP in S E100R0 was almost twice that in S E0R100 (Fig. 3b, NEP in Table 4).This is attributable to the difference in RH through RE, even though RH in S E100R0 was only 13% larger than that in S E0R100 (RH in Table 4).Cumulative NEP strongly affected ECP CNEP .The ECP CNEP of S E0R100 and that of S E100R0 were the longest (23 years) and the shortest (10 years), respectively, showing the same order as cumulative NEP (Fig. 3b, ECP CNEP in Table 4).However, ECP CNECB showed the opposite order: ECP CNECB of S E100R0 and S E0R100 were the longest (34 years) and shortest (23 years), respectively (Fig. 3c, ECP CNECB in Table 4).

Discussion
According to our VISIT model simulations of cool temperate forests, the effect of disturbance such as clear-cutting, conversion, and thinning on the long-term trend of NEP was greater than that of climate variation, even several decades after clear-cutting.Ito (2012) also reported that the effect of disturbance on NEP is larger than that of climate factors in deciduous broad-leaf forests in Japan.We found that interannual variation in the carbon balance was primarily attributable to climate variation.Thus disturbance controlled the long-term trend of the carbon balance, whereas climatic factors controlled the yearly variation.
The relationship of the young larch forest monthly NEP between model estimation and observation was not strong (Table S2 in the Supplement), and the simulated annual NEP of young larch forest was slightly biased compared with the observed NEP (Table S2 in the Supplement).One possible reason for these trends is that there might be differences in some ecological parameters between young and middleaged larch forests, whereas we used the same parameters for stands of both ages.The weak relationship and bias might have caused larger uncertainty in the case of the young larch forest than with the mature larch forest.Therefore, in the scenario test (Fig. 2, Table 3), the effects on carbon fluxes in the initial stage of recovery (i.e., in the first decade after disturbance) might have had some uncertainty.
The correlations of daily and monthly NEP were lower than those of GPP and RE (Table S2 in the Supplement), as reported in previous studies (Ichii et al., 2010(Ichii et al., , 2013;;Ueyama et al., 2010).NEP represents the small difference between large gross fluxes such as GPP and RE.Therefore, even a small variation in GPP and RE can have a large effect on the resultant NEP, causing low correlation of NEP.The low correlation of monthly NEP (R 2 = 0.42) and small slope for monthly RE (0.77) of the mature mixed forest might be attributable to the limitations imposed by a 1-year data set.Therefore, in the scenario test (Fig. 2, Table 3), there might have been larger uncertainty in NEP than in GPP or RE, and the carbon fluxes of S non−conv might have included larger uncertainty than in other scenarios.

The role of disturbance
In most cases, observation using the eddy covariance method detects carbon fluxes at one point within a long temporal sequence.Although the VISIT model could successfully simulate the carbon balance during several decades after disturbance without correct information about the amount of residues in the initial condition (Fig. 3a), residual biomass strongly affected the carbon balance in the first decade after disturbance (Fig. 3, Table 4).When we evaluate carbon debt and compensation periods for proper carbon management, we must consider not only the short-term carbon balance during the observation period but also the long-term sequential carbon balance after disturbance.
The carbon balance in the initial period, which is strongly affected by the amount of residues, is important for long-term sequential carbon fluxes.In addition, even several decades after the disturbance, the effect of clear-cutting on NEP was ascribed to decreased RH because the production of dead biomass decreased; this was consistent with the findings of Noormets et al. (2012).Therefore, treatment of residues plays an important role in long-term carbon balance and carbon management.Using a Monte Carlo simulation in a forest ecosystem in Washington State, Janisch and Harmon (2002) showed that the length of the NEP recovery period, which ranged from 0 to 57 years after disturbance, depended on the amount of residues.By applying the Biome-BGC model to a forest damaged by harvest and wildfire, Thornton et al. (2002) also showed that the NEP recovery period and carbon lost during this period depended on the amount of woody debris remaining on the site after the disturbance.Whereas Janisch and Harmon (2002) and Thornton et al. (2002) simply examined the effect of disturbance without taking climate variations into account, we simulated carbon balances while accounting for both disturbance and climate variations.Despite the difference in simulation design and the type of model, the results of all three studies suggest the importance of residues to the carbon balance of forest ecosystems.
Thinning is useful for enhancing NEP by decreasing RH, especially in middle-aged forests.By applying a carbon budget model to eastern Canadian red spruce forests, Taylor et al. (2008) found that thinning enhanced carbon assimilation because of increased productivity and reduced decomposition rate; this is consistent with our results.Dore et al. (2012) showed that thinning stimulated carbon assimilation of a www.biogeosciences.net/11/5139/2014/Biogeosciences, 11, 5139-5154, 2014 ponderosa pine forest in northern Arizona 3 or 4 years afterward, because it mitigated drought stress and allowed for carbon uptake.They also noted that it was difficult to detect the effect of thinning by using eddy covariance measurements; this is similar to the observational results of Scots pine forest presented by Vesala et al. (2005).According to Vesala et al. (2005), the reason why CO 2 flux was unaltered by thinning was that increased carbon absorption canceled out increased RH; this was opposite to our results.In their study, the amount of logging waste on the forest floor from thinning increased, whereas the effect of additional litter produced by thinning was small in our simulation.This could be one of the reasons why our results differed, and it suggests that treatment of residues affects carbon flux after thinning.Although our simulation clearly showed the effect of thinning, it might be difficult to detect by using field observations.

The role of climate factors
Variations of temperature and precipitation especially affected the interannual variation of annual NEP.Both factors directly affect GPP (or NPP) and RH because of their effects on photosynthesis and respiration, whereas solar radiation and VPD directly affect only NPP.However, RH is indirectly influenced by solar radiation and VPD by biomass change, because RH is influenced by the amount of dead biomass, which is influenced by NPP via the change in live biomass.
In the late 2000s, NEP decreased as temperature increased because enhanced RE would have been larger than enhanced GPP.Ito (2012) reported that, in a deciduous broad-leaf forest in central Japan, temperature was the meteorological factor with the largest effect on NEP; this was similar to our results.However, in that study, increasing temperature increased NEP, whereas we found the opposite.In both studies, both GPP and RE increase with increase in temperature.
In the case of our study, because increase in GPP exceeds increased in RE, NEP increases, and vice versa in the case of Ito (2012).Increased temperature affects many ecosystem processes (e.g., photosynthesis, plant respiration, soil respiration, biomass, phenology), and the sensitivities of all of these differ among plant species, soil types, and climate zones (Luo, 2007).Therefore, it is quite difficult to find a consistent pattern in the response to increased temperature (Luo, 2007), and further case studies should be conducted.
In our study, reducing the annual precipitation critically decreased annual NEP, as suggested by the scenario experiments.In contrast, Ito (2012) reported that precipitation did not affect NEP in a deciduous broad-leaf forest in central Japan.In the Ito (2012) study, the minimum annual precipitation was 1065 mm in 1994 (study period 1999-2009), whereas that in our study was 694 mm in 1984 (study period 1948-2010).Therefore, Ito (2012) did not detect severe drought.
A sensitivity analysis using the Biome-BGC model and 3 years of data (2001)(2002)(2003) at the Tomakomai site have shown that spring temperature and summer solar radiation are the primary factors affecting NEP (Ueyama et al., 2010), as also noted in our study, but precipitation did not influence NEP in the study by the Ueyama group.This is likely because Ueyama et al. (2010) added biases for each season, whereas our data included a year of severe drought in 1984 (694 mm), which had the lowest precipitation in 6 decades.This is likely to be the reason why we, but not Ueyama et al. (2010), noted a critical effect of precipitation on carbon balance.
In monsoonal Asia, it might be difficult to detect reduced carbon uptake caused by severe droughts on the basis of observations, because there are fewer severe drought events (the long-term mean annual precipitation from 1948 to 2010 was about 1000 mm at the Tomakomai and Teshio sites) than in North America or Europe.Saigusa et al. (2005) noted that moderate drought conditions may not reduce GPP but enhance it because of the prolonged growing season.However, our results suggested that severe droughts strongly affected the carbon balance, similar to the case in tropical forests in East Asia (Saigusa et al., 2008), and the effect of the severe drought in 1984 was larger than that of thinning.Although it is well known that drought affects soil respiration (Borken et al., 2006;Kopittke et al., 2013;Wang et al., 2014), reduced RH was not detected by our simulation under severe drought conditions.At the Tomakomai site, the relationship between soil respiration and soil water on the basis of chamber measurements was not clear (Liang et al., 2010).A much longer observation period is needed for detecting signals when severe drought occurs.Data based on a longer observation period are also useful for validating the VISIT model.
By applying an ecosystem model to future climate scenarios in aspen and jack pine forests in Saskatchewan, Canada, Wang et al. (2012) projected higher NEP under wet conditions than under dry conditions; this is consistent with our results.Yi et al. (2013) analyzed the effect of climate and fire on the carbon flux in a boreal and arctic ecosystem from 2000 to 2010 by using field and satellite observations and a satellite-oriented model.They found that NEP recovery after wildfire was enhanced under warmer conditions and reduced under drought conditions; this is similar to our results.However, they also reported that the effects of drought and temperature on NEP were larger than that of wildfire in most boreal and arctic areas; this was the opposite of our results.Yi et al. (2013) focused on the regional scale (boreal and arctic area), whereas we focused on the point scale in a cool temperate forest.The effect of disturbance depends on the area and strength of disturbance because various types of disturbance occurred heterogeneously at the regional scale.The effect of climate might also differ among regions.In future research, we need to evaluate the spatial effect of disturbance and climate by gathering information on disturbance in Hokkaido or Japan.

Ecosystem compensation points
The simulated ECP NEP at the Tomakomai site (5 years) was close to the observed ECP NEP calculated by using eddy flux data.ECP NEP was 7 years at the Teshio site.ECP CNEP and ECP CNECB were 8-34 and 20-91 years, respectively, on the basis of observational and previous data for the Teshio site (Aguilos et al., 2014).The simulated ECP CNEP (14 years) and ECP CNECB (23 years) at the Tomakomai site were within these ranges.According to Aguilos et al. (2014), ECP NEP in disturbed forests in cool temperate or boreal regions ranged from 7 years (clear-cut boreal upland forest in Canada; Howard et al., 2004) to 20 years (wildfire-damaged boreal forests in North America; Amiro et al., 2010;Grant et al., 2010), and ECP CNEP ranged from 3-17 years (clear-cut boreal upland forest in Canada; Howard et al., 2004) to 11-92 years (wildfire-damaged boreal forest in Canada; Goulden et al., 2011).Both the ECP NEP and the ECP CNEP values at the Tomakomai and Teshio sites were smaller than these values or at the lower limit, because larch forests have a high photosynthetic capacity (Hirata et al., 2007).The maximum GPP at light saturation at the TMK site was about 45 µmol m −2 s −1 ; this was higher than in a larch forest in Siberia (about 5 µmol m −2 s −1 ), an evergreen needle-leaf forest in Japan (about 20 µmol m −2 s −1 ), and an evergreen broad-leaf forest in Thailand (about 35 µmol m −2 s −1 ) (Hirata et al., 2008).Another reason is that the site's dense understories contributed to a large GPP in the initial period after disturbance (Takagi et al., 2009;Aguilos et al., 2014).According to the long-term simulation of the Tomakomai site (Fig. 2), the GPP of the understory accounted for 75 % of the total GPP in the first decade after disturbance (13 and 10 t C ha −1 year −1 for total GPP and understory GPP, respectively).Consequently, the recovery period tended to be shorter than in other forest types.
The ecosystem compensation point is controlled by the amount of residues and annual NEP after disturbance.As residual biomass increases, the ecosystem compensation point becomes longer.As annual NEP increases, the ecosystem compensation point becomes shorter.Therefore, not only ecosystem productivity but also residues related to forest practice are key factors in the carbon recovery period.

Conclusions
We evaluated the long-term effects of climate and disturbances on the carbon balance of forest ecosystems in a cool temperate region by using the process-based model VISIT.Clear-cutting strongly affected the carbon balance, not only just after clear-cutting but also 52 years after the disturbance.The effect of clear-cutting was larger than that of climate, even 52 years later.Disturbance controlled the long-term trend in carbon balance, whereas climate factors controlled the yearly variation.Among the climate factors, increased temperature and severe drought played vital roles in interannual variation of the carbon balance.The three ecosystem compensation points are useful indices for evaluating the recovery of a forest ecosystem from disturbance with respect to the carbon balance.Our findings suggest that the carbon balance of a forest ecosystem is controlled by the amount of residues, such that the recovery period can be controlled by the management of residues.The amount of residues strongly affected carbon release just after clear-cutting as well as the recovery period for NEP.W scalar is the minimum of two suppression terms: soil moisture and soil air space.Details of the functions are given by Ito and Oikawa (2002).
In the case of deciduous vegetation, such as larch forests, phenology is an important process.In the scheme of VISIT, there are four phenological periods: dormancy (from defoliation to leaf-out), emergence of new leaves, vegetative growth, and defoliation.The number of days between the beginning and end of leaf flush is determined by the growing degreedays, which is the cumulative temperature above 1 • C. The number of days between the beginning and end of leaf shedding is determined by cumulative coolness (degree-days below 12 • C).In the case of evergreen vegetation, the entire period is the vegetative growth period.Details of the phenology and turnover are given by Ito and Oikawa (2002) and Ito et al. (2005).
The process of litter fall or turnover brings carbon in the live vegetation compartments to the soil litter compartments.The litter fall (L) of each vegetation compartment is represented as follows: where SL is the specific litter fall rate or turnover rate per unit biomass, except for seasonal leaf abandonment of deciduous forest simulated in the phenology scheme.
Daily NPP is the difference between GPP and RA, whereas daily NEP is the difference between NPP and RH.In our study, positive NEP represented carbon uptake by the ecosystem and negative NEP represented carbon release from the ecosystem.The biophysical process of carbon allocation, hydrological processes, and meteorological calculations are described by Ito and Oikawa (2002).
When we simulated a disturbance such as clear-cutting or thinning, all or some of the carbon in the vegetation compartments was exported out of the ecosystem and the residues were moved to the litter compartments.When we take exported carbon into account, it is necessary to use the NECB, which is NEP minus the exported carbon, so we performed this calculation.Because clear-cutting is often accompanied by land-use conversion (Adachi et al., 2011), ecophysiological parameters were changed from those of the former vegetation type to those of the subsequent one when a conversion event occurred.Representative ecophysiological parameters of mixed and deciduous needle-leaf forests and understory are presented in Table A1.
The VISIT model has been validated by using several flux data sets covering tropical to subarctic biomes, including mature mixed forest in Teshio and middle-aged larch forest in Tomakomai (Ito, 2008;Ichii et al., 2010Ichii et al., , 2013)).In these studies, the parameters in Table A1 were determined manually by using a trial-and-error method to fit the output values of NEP, GPP, RE, NPP, and biomass to those derived from observation.In this study, only P max * of deciduous needleleaf forest was modified with reference to the maximum GPP at Tomakomai, as obtained by using the eddy covariance method (Hirata et al., 2007).
The Supplement related to this article is available online at doi:10.5194/bg-11-5139-2014-supplement.

Figure 1 .
Figure 1.Seasonal and interannual variations in net ecosystem production (NEP), gross primary production (GPP), and ecosystem respiration (RE) at the Teshio (TSE) and Tomakomai (TMK) sites.Black lines, output of the VISIT model; grey lines, observed values.A 2-week moving window was applied to the daily values of the model and observations. Figure2

Figure 2 .
Figure 2. Effects of climate, CO 2 , thinning, land-use change, and disturbance on (a) annual net ecosystem production (NEP), (b) NEP ( NEP), (c) gross primary production ( GPP), (d) ecosystem respiration ( RE), and (e) heterotrophic respiration ( RH) at the Tomakomai site.The values for climate, CO 2 , clearcutting, and conversion were calculated by subtracting each flux of S const−climate , S const−CO2 , S non−cut , and S non−conversion from each flux of S non−thin .The values for thinning were evaluated by subtracting each flux of S const−thin from each flux in the full scenario (S full ).Positive values represent the enhancement effects on each flux, and negative values indicate the opposite.
by strong winds in typhoon Toyamaru in 1954.Thinning was  conducted in 1970Thinning was  conducted in  , 1985Thinning was  conducted in  , and 2004.   .

Table 2 .
Average annual carbon fluxes estimated for three forest types -mature mixed forest (2002 at Teshio), young larch forest (2003-2012 at Teshio), and middle-aged larch forest (2001-2003 at Tomakomai) -by field observation and the VISIT model.Units are t C ha −1 year −1 .Values in parentheses are standard deviations.

Table 3 .
Average effects of climate, air temperature, precipitation, solar radiation, vapor-pressure deficit (VPD), CO 2 , clear-cutting, conversion, and thinning on net ecosystem production at the Tomakomai site for each decade after the disturbance in 1958.Positive values represent enhancement effects on each flux, and negative values indicate the opposite.Units are t C ha −1 year −1 .Values in parentheses are standard deviations.

Table 4
, which lists the annual values just after clear-cutting in 1958.Long-term effects are also shown in Table4, which gives the cumulative annual values of carbon fluxes from 1958 to 2010.

Table 4 .
Sensitivities of net ecosystem production (NEP), gross primary production (GPP), ecosystem respiration (RE), heterotrophic respiration (RH) (upper table) and ecosystem compensation point (ECP) (lower table) to changes in residues.In the upper table, the first row shows the values of NEP, GPP, RE, and RH immediately after clear-cutting in 1958.

Table A1 .
Representative parameters of the VISIT model for the mixed and deciduous needle-leaf forests and understory.