Interactive comment on “ Quantifying legacies of clearcut on carbon fluxes and biomass carbon stock in northern temperate forests ” by W .

1. This study employed a well-established process-based forest model (PnET-CN) to evaluate the effect of clearcut on carbon (C) flux trajectory in two widespread plant functional types (deciduous broadleaf forest and evergreen needleleaf forest) in the upper Midwest region of Wisconsin and Michigan. The trajectory analysis of C flux after clearcut makes this study quite interesting. Results suggested that harvest have a big influence on early stage of forest succession, but only had little effects on late stage. The method used in this study is solid, and results met with expected recovery

We compared the effects of stand-replacing harvesting on C fluxes and stocks using two chronosequences of eddy covariance flux sites for deciduous broadleaf forests (DBF) and evergreen needleleaf forests (ENF) in the Upper Midwest region of northern Wisconsin and Michigan, USA.The average values of normalized root mean square error (NRMSE) and the Willmott index of agreement (d ) between simulated and inferred from observation variables including gross primary productivity (GPP), ecosystem respiration (ER), net ecosystem productivity (NEP), LAI, and AGC in the two chronosequences were 20 % and 0.90, respectively.Simulated GPP increased with stand age, reaching a maximum (∼ 1200-1500 g C m −2 yr −1 ) at 11-30 years of age, and leveled off thereafter (∼ 900-1000 g C m −2 yr −1 ).Simulated ER for both forest types was initially as high as ∼ 700-1000 g C m −2 yr −1 in the first or second year after clearcuts, decreased with age (∼ 400-800 g C m −2 yr −1 ) before canopy closure at 10-25 years of age, and increased to ∼ 800-900 g C m −2 yr −1 with stand development after canopy recovery.Simulated NEP for both forest types was initially negative with the net C losses of ∼ 400-700 g C m −2 yr −1 for 6-17 years after harvesting, reached the peak values of ∼ 400-600 g C m −2 yr −1 at 14-29 years of age, and became stable and a weak C sink (∼ 100-200 g C m −2 yr −1 ) in mature forests (> 60 years old).The decline of NEP with age was caused by the relative flatting of GPP and gradual increasing of ER.ENF recovered slower from net C source to net sink and lost more C than DBF, suggesting ENF are likely slower to recover C assimilation capacity after stand-replacing harvests due to slower development of photosynthesis with stand age.Model results indicated

Introduction
Disturbance has been widely recognized as a key factor influencing ecosystem structure and function at decadal to century scales (Magnani et al., 2007;Williams et al., 2012;Kasischke et al., 2013).Harvest is an important anthropogenic disturbance shaping North American forest landscapes.Approximately 61 000 km 2 of forests were affected by harvests every year during the 2000s (Masek et al., 2011).Harvests affect forest age structure and alter the forest carbon (C) balance (Magnani et al., 2007;Pan et al., 2011;Williams et al., 2012).Quantifying the legacy of harvest disturbances is essential for predicting forest C dynamics, informing climate policy-making, and improving forest management under the context of climate change.Here, we focus on assessment of an ecosystem model C cycle response from one type of harvest, the clear-cut.
Harvests transfer living biomass C to harvested wood C and litter C, resulting in successional changes in C fluxes and stocks.Leaf biomass increases rapidly in secondary succession and then typically stabilizes at a certain level that is determined by light, water, nutrient availability, and forest type (Sprugel, 1985).Gross primary productivity (GPP) thus increases gradually over time, reaches maximum in middle age, and in response to nutrient limitations and aging responses slightly declines thereafter (Odum, 1969;Chapin et al., 2002).The successional change in plant respiration (autotrophic respiration) after stand-replacing harvesting is similar to that of GPP, although C use efficiency (the ratio of net primary productivity to GPP, NPP/GPP) could decline with forest age (DeLucia et al., 2007).As a result, living tree biomass C gradually increases following a typical logistic growth curve (Odum, 1969;Sprugel, 1985) Heterotrophic respiration following stand-replacing harvesting could be stimulated at the beginning of stand development because the removal of trees alters the environmental conditions (e.g., soil temperature, moisture, and nutrients) and possibly leads to changes in litter quantity depending on harvesting types (e.g., stem-only harvesting) (Chapin et al., 2002).Heterotrophic respiration is expected to gradually decrease thereafter because the regrowing forest reduces light, water, and nutrient availability to the soil (Chapin et al., 2002) and the amount of decomposable soil organic matter from the prior forest and harvest residue (e.g., litter, coarse woody debris, and soil organic C) also gradually decreases.Over time, heterotrophic respiration could be enhanced again because of accumulation of woody debris and litter with stand development.This theorized successional trajectory in ecosystem respiration (ER; the sum of autotrophic and heterotrophic respiration) may also be strongly influence by harvest types and forest composition.Thus the trajectory of heterotrophic respiration (and consequently total ecosystem respiration) with age is not well-understood (Amiro et al., 2010).Observational studies to date have shown that forest ecosystems generally become C sources (negative net ecosystem productivity, NEP) immediately following stand-replacing harvests, approach the maximum NEP as they mature, and then experience a gradual decline in NEP thereafter (e.g., Law et al., 2003;Gough et al., 2007;Goulden et al., 2011), following the trajectories hypothesized by Odum (1969) and Chapin et al. (2002).
The changes of C fluxes and stocks after harvesting have been examined in many forest ecosystems using ecological measurements (e.g., eddy covariance observations) from chronosequences using a space-for-time substitution approach (e.g., Law et al., 2003;Gough et al., 2007;Goulden et al., 2011).The trajectories and amplitude of C fluxes and stocks vary with forest ecosystem types (Amiro et al., 2010).For example, Noormets et al., (2007) reported that a young red pine (Pinus resinosa) stand at 8 years of age was a net C sink (313 ± 14 g C m −2 yr −1 ), but a young hardwood site at age of 3 was a net C source (−128 ± 17 g C m −2 yr −1 ) over the growing season in northern Wisconsin, USA.Young stands in northern Wisconsin may become net C sinks within 10-15 years after harvesting (Noormets et al., 2007).More rapid recovery Introduction

Conclusions References
Tables Figures

Back Close
Full after stand-replacing harvesting (< 6 years) was found for temperate forests in northern Michigan (Gough et al., 2007).These studies have produced a wealth of information on ecosystem C dynamics after stand-replacing disturbances, and this information can be translated to more process-based and quantitative understanding of disturbance effects on the C cycle using ecosystem models (Goulden et al., 2011).Process models require evaluation on how source/sink transition and long-term carbon flux dynamics respond to differences in vegetation type, harvest intensity, and age since clearing.
Although using the chronosequence approach to evaluate the changes of ecological processes with age after disturbances is attractive, this approach is often limited by biological and climatic factors (Yanai et al., 2003;Bond-Lamberty et al., 2006) and the lack of full representation of stand development stages.Process-based ecosystem models provide a means of quantifying the effects of disturbances on C dynamics under changing climate over various spatial and temporal scales.Ecosystem models have been used to assess the effects of clearcuts and climate change on forest C dynamics at the stand/ecosystem (e.g., Bond-Lamberty et al., 2006;Grant et al., 2009;Wang et al., 2012b) or regional scales (Desai et al., 2007;Dangal et al., 2014).Moreover, ecosystem models can also be used to assess forest C dynamics under various scenarios of climate change and harvesting regimes (e.g., Albani et al., 2006;Peckham et al., 2012), since these models have been developed based on physiological, biogeochemical, and ecological theories.However, few studies have used ecosystem models to examine the changes of C fluxes and stocks with stand regrowth after standreplacing disturbances for forest chronosequences.
The objective of this study was to evaluate the ability of an ecosystem model to capture the trajectories of forest C dynamics after a stand-replacing disturbance in two northern forest chronosequences, examine by which processes successional trajectories evolve in these systems, and test the role of forest composition on successional trajectory of C fluxes.The two chronosequences, whose sites all were initiated by some sort of stand-replace harvest, can be broadly classified along two plant functional types: deciduous broadleaf forests (DBF) and evergreen needleleaf forests (ENF).We applied Introduction

Conclusions References
Tables Figures

Back Close
Full a process-based forest ecosystem model, PnET-CN (Aber et al., 1997;Ollinger et al., 2002), to simulate the effects of clearcut on forest C dynamics, and evaluated the simulated C fluxes and stocks for both forest types using in-situ measurements (e.g., eddy covariance observations and aboveground biomass C).We hypothesized that (1) both DBF and ENF will have similar successional patterns in C fluxes (GPP, ER, and NEP) and above-ground biomass C stocks (AGC) after stand-replacing harvests but (2) DBF will recover faster than ENF from a net C source to a net C sink and lose a smaller amount of C (negative NEP) following a stand-replacing harvest.

Study sites and data
Our study sites consist of 8 eddy covariance sites in the Upper Midwest region of northern Wisconsin and Michigan (Chen et al., 2008; Table 1).The study area has a humidcontinental climate, with hot summers and cold winters.The mean annual temperature is 4.4 • C and the mean annual precipitation is 768.9 mm (as measured between 1981 and 2010 at Rest Lake weather station, 46.12 • N 89.87 • W, http://www.ncdc.noaa.gov).
The dominant soil type is glacial sandy loam and loamy tills (Noormets et al., 2008).The region has been strongly influenced by forest industry.Most forest stands less than 100 years old in this region regenerated following harvesting operations (Amiro et al., 2010).
Our sites consist of four DBF sites (YHW, IHW, WIC, and UMBS) and four ENF sites (YRP, YJP, IRP, and MRP).The four DBF sites range from 3 to 86 years in age and constitute a chronosequence.The four ENF sites also represent a chronosequence with stand age ranging from 8 to 66 years.We obtained monthly C fluxes (observed NEP and its inferred data products GPP and ER) from AmeriFlux (http://public.ornl.gov/ameriflux/) for the eight eddy covariance flux tower sites (Table 1).Harmonized level 4 data were used in this study.These flux data have been described and used in Introduction

Conclusions References
Tables Figures

Model description
The PnET-CN model is a process-based forest ecosystem model that can simulate C, nitrogen, and water dynamics at monthly time steps.PnET-CN is driven by temperature, precipitation, photosynthetically active radiation (PAR), wet and dry nitrogen deposition, and atmospheric CO 2 concentration (Aber and Driscoll, 1997;Aber et al., 1997;Ollinger et al., 2002).The model has been applied and tested in the USA and Europe to simulate effects of climate variability, rising atmospheric CO 2 , ozone pollution, and disturbance on ecosystem processes and functions (e.g., Aber et al., 2002;Pan et al., 2009;Peters et al., 2013).One of the unique features of PnET-CN is to simulate potential photosynthesis using foliar nitrogen concentration and light use efficiency in a multilayered canopy (Aber and Federer, 1992).Actual GPP is constrained by air temperature, vapor presser deficit, and soil water availability.The effects of elevated CO 2 concentration on leaf photosynthetic rates are calculated using constant ratios of leaf internal to ambient CO 2 concentration (C i /C a ) (Ollinger et al., 2002).PnET-CN incorporates a total of seven C pools, five of which are structural C pools (foliage, woods, fine roots, woody debris, and soil organic matter) and two of which are non-structural C pools stored in woods and fine roots.Photosynthetic production is allocated to each living plant component (i.e., foliage, woods, and fine roots) and to growth and maintenance respiration.Living tree biomass is transferred to dead woody biomass and/or to soil organic C through leaf and root turnover, tree mortality, and disturbance.The decomposition of coarse woody debris is a constant fraction of its C content.The decomposition of soil organic C is calculated as a function of maximum decomposition rate and effects of temperature and soil moisture.Introduction

Conclusions References
Tables Figures

Back Close
Full PnET-CN includes a complete nitrogen cycle, and simulates nitrogen mineralization and nitrification, plant nitrogen uptake, allocation, and leaching losses.Nitrogen depositions are imposed into corresponding soil nitrogen pools (NH 4 and NO 3 ).As with C pools, nitrogen is divided into five structural pools (foliage, woods, fine roots, woody debris, and soil organic matter) and one non-structural nitrogen pool stored in the trees.
C and nitrogen cycles interact closely in the model.High leaf nitrogen concentration increases net photosynthesis rate -in the absence of water stress, resulting in the high demand for non-structural nitrogen in plant tissues (Aber et al., 1997).When plant non-structural nitrogen is low, plant nitrogen uptake efficiency from available soil mineral nitrogen is increased in the model (Aber et al., 1997).In addition, high C : N ratios in biomass, litter, and soil organic matter reduce net mineralization rates.In general, the nitrogen cycle in the model is governed by a negative feedback loop.
The model also simulates key hydrological processes including rainfall interception, evaporation, transpiration, surface runoff, and drainage at each time step.Rainfall interception is treated as a constant fraction of precipitation.Transpiration is estimated based on water use efficiency.Similarly, surface runoff is treated as a constant fraction of the difference between precipitation and evaporation.Drainage is estimated when potential soil water exceeds soil water holding capacity.
Prescribed disturbance events can be simulated in the model through four parameters: disturbance year, disturbance intensity, biomass removal fraction, and the loss rate of soil organic matter.In this study, when stand-replacing disturbance events occur, a uniform forest type was assumed to be regenerated on-site.Moreover, minimum leaf area index (LAI) of 0.1 was assumed to regulate maximum potential foliage mass that controls leaf production only in the first year after clearcuts.The photosynthetic production is transported to plant non-structural C pool where C could be allocated to leaves, stems, and roots.There is thus no need for initialization (e.g., stand density) after disturbances in the model.More details about the model structure and processes have been described elsewhere (Aber et al., 1997;Ollinger et al., 2002).Introduction

Conclusions References
Tables Figures

Back Close
Full

Model inputs
The model inputs include temperature, precipitation, PAR, wet and dry nitrogen deposition, atmospheric CO 2 concentrations, and disturbance history.The climate data used in all simulations were derived from the Daymet database (Thornton et al., 2012).
For each site, monthly maximum temperature, minimum air temperature, and precipitation were calculated from the daily Daymet data for the period 1981-2010.PAR (mol m −2 s −1 ) was estimated from solar radiation (RAD, MJ m −2 day −1 ) using the empirical relationship (PAR = 2.05 RAD) (Aber et al., 1996).Those measurements from 1981 to 2010 were repeated to create the time series from 1850 to 1980.Annual rates of wet and dry nitrogen deposition were obtained from the United States Environmental Protection Agency (EPA; http://java.epa.gov/castnet/clearsession.do).The nitrogen deposition rates were measured at the Wellston station (44.22 • N; 85.82 • W) for the period 1994-2011.We also obtained the nitrogen deposition rates in 1860 estimated by Galloway et al. (2004).For each year prior to 1994, we used an exponential ramp function to estimate the annual deposition rates by interpolating the historical (1860) and current nitrogen deposition rates.Monthly wet deposition rates, needed for the model, were generated from annual wet nitrogen deposition through the weighted coefficients (the ratio of monthly precipitation to total precipitation from March to November).We assumed that there is no wet nitrogen deposition in winter.The soil water holding capacity in the rooting zone (100 cm) for each site was derived from the gridded multi-layer soil characteristics dataset (STATSGO, Miller and White, 1998).
For each site, we prescribed the disturbance events using the site disturbance history (Table 1).For each stand-replacing harvest, stand mortality was assumed to be 100%.The merchantable wood removal (biomass removal out of the ecosystem) fraction was assumed to be 0.8 in this study.We also conducted a sensitivity analysis to these assumptions, as described below.The soil removal fraction was assumed to be zero, given that the content of soil organic C might not be considerably affected by harvesting (Johnson and Curtis, 2001;Yanai et al., 2003).Introduction

Conclusions References
Tables Figures

Back Close
Full

Parameterization, initialization, validation, and sensitivity analysis
PnET-CN has been parameterized and tested for temperate DBF (Aber et al., 1997;Ollinger et al., 2002;Peters et al., 2013), temperate ENF (Aber et al., 1997;Peters et al., 2013), and mixed forests (Aber et al., 1997) for forest productivity, net nitrogen mineralization, and foliar nitrogen concentrations.The parameter values used in this study are given in Table 2. To apply the model to the transient simulation period (1860-2010), a 200-year spin up run was conducted to ensure that the equilibrium (∆NEP < 10 g m −2 month −1 and ∆soil organic C < 1%) was reached for each chronosequence site.The climate normals , pre-industry nitrogen deposition rates, and historical CO 2 concentrations were used for the spin up runs.
To examine the stand-replacing harvest legacies, we conducted all simulations using the site disturbance history (Table 1), vegetation parameters (Supplement Table S1), climate, nitrogen deposition, and atmospheric CO 2 for each of the chronosequence sites.The model simulations were evaluated against C fluxes (GPP, ER, and NEP), AGC, and LAI data collected at the eddy covariance flux sites.We used two statistical measures to evaluate the overall model performance: the normalized root mean square error (NRMSE) and the Willmott index of agreement.The NRMSE (Eq. 1) was used to assess the difference between predicted (P ) and observed (O) variables.
where O max and O min are the maximum and minimum observed values, respectively; i is the i th observation; and n is the total number of observations.A value close to 0 as: (2) A value of 1 indicates a perfect agreement and a value near 0 indicates a weak agreement (Willmott, 1982).
Model sensitivity to changes in harvesting practices on ecosystem C dynamics during the secondary succession was conducted.The model was run at WIC and MRP for 100 years after scenario harvests in 1910, using the same climate data sequence.Sensitivity scenarios involved applying the stand mortality (80% and 60%, compare to 100% in the model test) to reveal effects of different harvest intensity scenarios on C dynamics.

Evaluation of modeled carbon fluxes and stocks
The simulated C fluxes were generally consistent with eddy covariance derived C fluxes for both DBF and ENF sites (Figs. 1 and 2).The NRMSE between simulated and tower fluxes (GPP, ER, and NEP) were between 10-21% (Table 2) Overall, the simulated and observed stand characteristics (LAI and AGC) showed good agreement (Table 2 and Fig. 3).The model slightly underestimated LAI for the young forest sites, and overestimated LAI for the mature forest sites.Generally, the model overestimated AGC for the mature forest sites, compared to the observations.The NRMSE of simulated and observed AGC was 28% and 31% for LAI.The Willmott index of agreement between simulated and observed AGC and LAI was 0.95 and 0.97, respectively.

Legacy of clearcut on carbon fluxes and stocks
PnET-CN generally captured the changes of C fluxes following the clearcuts for each chronosequence site (Fig. 4).The predicted annual GPP generally increased with time since disturbance and approached the peak values (∼ 1200-1500 g C m −2 yr −1 ) between 11 and 26 years of age and between 29 and 30 years of age for the DBF (IHW, WIC, and UMBS) and the ENF (IRP and MRP) sites, respectively; thereafter, the forest stands reached maturity and GPP became relatively stable with mean values of 940-1000 g C m −2 yr −1 .
Predicted annual ER was initially as high as 860-1030 and 710-860 g C m −2 yr −1 within the first two years for the DBF and the ENF sites, respectively.During canopy recovery, predicted ER generally decreased to 620-780 g C m −2 yr −1 between 10 and 25 years of age for the DBF sites and to 360-380 g C m −2 yr −1 between 14 and 17 years of age for the ENF sites (Fig. 4).For forest ages older than 60 years, the predicted annual ER for both forest types showed a relatively flat pattern, contrary to theoretical expectations, arising from the little change of both autotrophic and heterotrophic respiration with age (Supplement Fig. S1).Average annual ER for mature forests was 810-880 and 780 g C m −2 yr −1 for the DBF sites (WIC and UMBS) and the ENF (MRP) site, respectively.
As expected, the ratio of annual GPP to annual ER (GPP : ER) simulated by PnET-CN was low during the early years after clearcutting for both forest types (Fig. 5).Within ∼ 6 years for the DBF sites and ∼ 17 years for the ENF sites, the GPP : ER ratio grad-8800 Introduction

Conclusions References
Tables Figures

Back Close
Full ually increased and its average value became larger than 1 (NEP > 0).The simulated peak GPP : ER ratio for DBF (1.6) occurred at 18 years after stand-replacing harvests, and the simulated peak ratio for ENF was 1.8 at 26 years.After those peaks, the ratio became relatively stable, with the mean values of 1.1 and 1.2 for mature DBF and mature ENF, respectively.
The model predicted negative NEP (C source) for the 6 and 17 years after standreplacing harvests for the DBF and the ENF, respectively (Fig. 4).The simulated peak annual net C loss occurred in the first or second year after clearcutting.The average C loss was 530-710 g C m −2 yr −1 for the DBF sites and 380-400 g C m −2 yr −1 for the ENF sites.The total C loss was 3.2-4.3and 6.4-6.9 kg C m −2 for the DBF and the ENF sites, respectively.The maximum net C gain was 387-433 g C m −2 yr −1 at 14-26 years of age for the DBF sites (WIC and UMBS) and was 567-602 g C m −2 yr −1 at 29 years of age for the ENF sites (IRP and MRP).Simulated annual NEP decreased thereafter and became as low as 120-180 g C m −2 yr −1 after 17-31 years for the DBF sites and 170 g C m −2 yr −1 after 44 years for the ENF sites.
Forest canopy as measured by LAI gradually recovered over time following clearcuts (Fig. 6).LAI fully recovered within 10-15 years after disturbance for the DBF sites and within 40 years of age for the ENF sites.The gradual recovery of LAI led to the gradual increase in GPP and the subsequent increase in AGC (Fig. 7).In general, AGC recovered much more slowly than C fluxes and LAI.The changes of simulated AGC followed the logistic growth curve with slow accumulation in the early years, fast accumulation in the middle age, and slow accumulation afterwards.The predicted LAI and AGC generally fell within the range of observed values across two chronosequences (Figs. 3, 6,  7).For mature forests (> 60 years of age) in 2010, the DBF sites generally stored more C in aboveground biomass than the ENF sites (10-12 vs. 8.5 kg C m −2 ; Fig. 7).

Sensitivity to harvest
Harvest intensity had little effect in the long term C dynamics for both forest types, but has obvious effects during early succession (Fig. 8).Greater harvest intensity led

BGD Introduction Conclusions References
Tables Figures

Back Close
Full to earlier rising GPP (Fig. 8a and f) and LAI (Fig. 8d and i), but delayed reduction in ER (Fig. 8b and g), finally resulting in latter rising NEP (Fig. 8c and g).High harvest intensity (e.g., 100% removal of living trees) also directly reduced living tree AGC (Fig. 8e and i).By reducing harvest intensity parameter to 80% and 60% from 100% used in the original model, average annual NEP over 100 years for DBF decreased by 104 and 88 g C m −2 yr −1 , respectively.The increased remaining tree biomass resulted in increasing AGC about 12% and 16%, respectively after a 100-yr harvest cycle.For ENF average annual NEP decreased about 1% and AGC decreased nearly 6% for both reduced harvest intensity scenarios.

Carbon fluxes and stocks following stand-replacing harvests
Overall, the model generally simulates C flux post-harvest trajectories that would be expected for LAI, NEP, GPP, autotrophic and heterotrophic respirations.The model showed that after forest recovery, expected increases in ER with age do not occur, in according with observations at the chronosequence sites.However, the model was unable to simulate high GPP rates estimated by eddy covariance in mature forests, regardless of vegetation type, suggesting potential improvements for model simulation of secondary succession, as discussed below.
Our simulations show that LAI increased rapidly first and then stabilized during the following development stages, given that the model estimates foliage growth through the parameter of maximum relative growth rate (Table 2) with the limitation of current foliage biomass and resource availability.This modeled response is consistent with earlier observations that foliage biomass increases rapidly after disturbance and then stabilizes (Sprugel, 1985).Our chronosequence-based results are generally consistent with previous results.For example, Goulden et al. (2011) reported that LAI along a chronosequence of boreal forest stands increased rapidly from 0.3 m 2 m −2 1 year Introduction

Conclusions References
Tables Figures

Back Close
Full The simulated successional change in annual GPP for both forest types generally follows the trajectory hypothesized by Odum (1969).However, while a slight decrease in GPP is hypothesized in Odum's trajectories, our simulations show a relatively flat GPP in mature forests (Figs. 4 and 10).In the model, GPP tracks LAI in the absence of significant changes in light, water or nutrient stress.As LAI stabilizes in mature forests, GPP does too.Our results are consistent with previous studies showing relatively flat pattern in GPP with age after 20 years following harvests in temperate pine forests in Florida (Clark et al., 2004), northern temperate DBF in Wisconsin (Desai et al., 2008), and boreal pine forests in Saskatchewan (Zha et al., 2009).Furthermore, Humphreys et al. (2006) reported continuous increases of GPP with increasing forest age for coniferous forests using three different stands at different stages of development (2, 14, and 53 years of age) following clearcuts.However, northern temperate ENF shows a small difference in GPP between young and mature sites (Noormets et al., 2007;Desai et al., 2008).Site-to-site variations in species and soil fertility could result in the differences in GPP successional trajectory after clearcuts that deviate from expected theory or model mechanisms.Also, the C flux and stock trajectories derived from chronosequence studies (both ecological measurements and eddy covariance observations) could be limited to some extent given the fact that the chronosequence method could not encompass a full range of forest development stages (e.g., Clark et al., 2004;Humphreys et al., 2006;Noormets et al., 2007).We lack significant oldgrowth data, though Desai et al. (2005) show that a nearby old-growth mixed forest did have slightly lower GPP and significantly higher ER than nearby DBF.
Our simulations show that annual ER for secondary temperate forests slightly declines in the first ten years because autotrophic respiration is low at first after the removal of trees.Amiro et al. ( 2010) reported that ER reduced in the very first year Introduction

Conclusions References
Tables Figures

Back Close
Full following harvests for a number of eddy covariance flux sites over North America.
Previous field studies showed that ER following clearcuts increased with forest age (e.g., Humphreys et al., 2006;Zha et al., 2009), partly supporting our result that ER slightly increases after the short decline period (10-25 years of age) in northern temperate forests until they reach maturity.Martin and Bolstad (2005) reported that chamber-based soil respiration in DBF of northern Wisconsin ranged from 857-1143 g C m −2 yr −1 in 1998 and 1013-1357 g C m −2 yr −1 in 1999, which is higher than tower ER (825 ± 133 g C m −2 yr −1 , WIC) from 1999-2006 in the same region.Soil respiration of 690 g C m −2 over the growing season of 2005 in a mature DBF near WIC tower site was reported (Tang et al., 2009).Our simulated respiration components (e.g., soil respiration) for DBF were lower than those reported values (Supplement Fig. S1).The model underestimated GPP but well estimated ER for mature DBF sites compared to derived fluxes from tower, indicating the model might underestimate root autotrophic respiration.Eddy covariance derived ER were usually lower than Chamberbased estimation in the WIC site due to uncertainties induced by measurement methods, decoupling of surface and canopy fluxes at night, and spatial scaling (Bolstad et al., 2004;Cook et al., 2004).For ENF the model underestimated ER compared to tower derived ER (340 ± 96 g C m −2 yr −1 ).This could result from overestimated soil decomposition rate in the model.Our simulations also show that DBF has slightly higher soil respirations than ENF (Supplement Fig. S1), which is consistent with the finding that chamber-based soil respiration was slightly higher for DBF than for ENF (no significant difference) during 1999 and 2000 growing season in Wisconsin (Euskirchen et al., 2003).Variability in changes of ER in secondary forests after clearcutting differ because of the quantity and quality of soil organic C and litter C as a result of site conditions and harvesting types (e.g., Tang et al., 2009).
The trajectory of our simulated GPP : ER ratio is similar to the curve derived by Amiro et al. ( 2010) using EC observations and forest age from chronosequences across North America (GPP: ER = 1.23 • [1 − exp(−0.224• AGE)]).Our simulated ratios are within the observed range of 0.9-1.6 for the DBF sites (Fig. 5a), although the model might un-Introduction

Conclusions References
Tables Figures

Back Close
Full derestimate the ratios for mature sites according to derived GPP and ER for the two chronosequence sites.The growing season GPP : ER ratios are normally higher than the annual ratios because winter soil organic C decomposition is important to annual C balance (Aanderud et al., 2013).However, the simulated ratios for the ENF sites are much lower than tower-derived growth season ratios (1.9-4.7,Fig. 5b), and close to the annual range of 1.6-2.2estimated by Desai et al. (2008).The standard gap-filling methods may lead to overestimation of net ecosystem change due to the lack of winter C flux data for the ENF sites and two of the DBF sites (YHW and IHW).
Our simulated successional dynamics of NEP following clearcuts generally supported the trajectory hypothesized by Odum (1969) and Chapin et al. (2002).The similar trajectories, however, are induced by different causes.Odum's trajectories show declining GPP and relatively flat community respiration with time.Our simulated decline in NEP results from relatively flat GPP and growing ER with stand development (Figs. 4 and 8).This has been observed for northern temperate hardwood chronosequence sites (Desai et al., 2008), northern temperate pine forests (Peichl et al., 2010), and boreal DBF forests (e.g., Goulden et al., 2011).The recent North American Carbon Program (NACP) synthesis study showed a similar change in NEP based on eddy covariance chronosequence measurements across North America (Amiro et al., 2010).NEP can also be calculated as the difference between NPP and heterotrophic respiration.Chapin et al. (2002) hypothesized that heterotrophic respiration is initially high, decline in middle succession, and rise thereafter, while NPP reaches a peak in middle age and declines in old stands.The simulated successional trajectories in heterotrophic respiration support the pattern hypothesized by Chapin et al. (2002), whereas our simulated NPP does not decline in mature stands (Supplement Fig. S1).Previous studies support our simulated trajectory in heterotrophic respiration.For example, Pregitzer and Euskirchen (2004) reported that heterotrophic respiration was high (mean value of 970 g C m −2 yr −1 ) in young temperate forests, declined with age in middle succession, and increased with time for mature forests, although for old temperate forests (> 120 years) the decline in NPP reduced heterotrophic respiration.The decline of NPP Introduction

Conclusions References
Tables Figures

Back Close
Full with age is not predicted in this study.Validation of the NPP was not possible in this study due to the lack of NPP data available across all sites.Our simulated heterotrophic respiration for mature DBF is close to the observation of 502 ± 86 g C m −2 yr −1 in a mature DBF near UMBS tower site between 1999 and 2003 (Gough et al., 2008).However, ENF chronosequence sites in this study show NEP continually increasing with age.It is induced by relatively flat and low ER (340 ± 96 g C m −2 yr −1 ) and increasing GPP.
Although our model underestimated NEP and GPP for both the DBF and ENF sites in the Upper Midwest region (Figs. 1 and 2), our predicted NEP was comparable to estimates from other studies in a similar region.For example, our predicted maximum NEP for the ENF sites (567-602 g C m −2 yr −1 , 29 years of age) was slightly lower than the estimates (690 g C m −2 yr −1 , 15-20 years of age) for afforested white pine (Pinus strobus) forests in Ontario (Coursolle et al., 2012).For a northern temperate forest chronosequence study in northern Michigan, NEP higher than 200 g C m −2 yr −1 in young DBF forests could be derived from the reference forest (153 ± 115 g C m −2 yr −1 , 85 years of age) (Gough et al., 2007), suggesting our predictions (390-430 g C m −2 yr −1 , 14-26 years of age) for the DBF sites could be in the reasonable range.Furthermore, our predicted mean annual NEP (123-177 g C m −2 yr −1 ) for mature DBF sites (> 60 years) was close to estimates for other northern DBF, including a northern hardwood forest of central Massachusetts (200 ± 40 g C m −2 yr −1 , Barford et al., 2001), and four eastern North America DBF (167-236 g C m −2 yr −1 Curtis et al., 2002).
We found that the simulated AGC during forest regrowth gradually increased following the typical logistic growth curve (Sprugel, 1985).In the model, low NPP in the early stages results in slow AGC accumulation.Once the amount of NPP approximately equals annual dead biomass C that is largely controlled by the wood turnover rate, the trajectory of AGC reaches a plateau.Previous chronosequence studies have reported that AGC increased with increasing age (e.g., Peichl and Arain, 2006;Goulden et al., 2011;Powers et al., 2012).Powers et al. (2012) reported that AGC increased rapidly with age in young red pine stands across a chronosequence in northern Minnesota, USA.Very few young stands could limit their finding (Powers et al., 2012).Introduction

Conclusions References
Tables Figures

Back Close
Full

Differences between DBF and ENF
Our simulations showed clear ecosystem differences between the DBF and ENF sites, suggesting that DBF may reach a peak in LAI and GPP faster than ENF after clearcutting.More rapid recovery of LAI and GPP for DBF sites lead to sooner restoration for NEP and AGC, supporting our second hypothesis.The foliage related parameters such as FolRelGroMax and AmaxB mainly govern the differences in successional trajectories between the two forest types (Table 2).DBF is assumed to have more productive foliage than ENF, and more photosynthetic production then could lead to more foliage production.With this positive feedback in the model, GPP, NEP, and AGC of the DBF sites recover more rapidly than those of the ENF sites.Our findings are consistent with the chronosequence study reporting that the temperate DBF in northern Michigan rapidly became a net C sink after six years following disturbances (Gough et al., 2007), and ENF stands in northern Wisconsin became net C sinks within 10-15 years after harvesting (Noormets et al., 2007).Through analyzing United States Forest Service Forest Inventory and Analysis data, Williams et al. (2012) suggested that faster growth in AGC at high productivity sites caused higher C fluxes and stocks.Our findings are also consistent with a recent modeling study suggesting that temperate DBF switches positive NEP faster than temperate ENF after clearcuts, and DBF has a higher peak in NEP compared to ENF (Peckham et al., 2012).A modeling study conducted in boreal forests also reported that low productive boreal ENF needed 1-3 years more to attain a positive NEP than boreal DBF after clearcuts in Saskatchewan, Canada (Wang et al., 2012b).These observed and modeled successional changes further reflect that DBF tend to have higher photosynthetic capacity compared to ENF in the early stage of stand development following stand-replacing harvests.

Model limitations and challenges
PnET-CN can explicitly simulate the effects of disturbance, pollution, and climate change conditions on forest C dynamics (e.g., Ollinger et al., 2002;Pan et al., 2009; Pe-

Conclusions References
Tables Figures

Back Close
Full  , 2013).Despite the capability of the model, we do recognize that the model has some limitations in simulating harvesting effects, and accurate representation of the trajectories of C fluxes and stocks following harvests still remains a challenge.First, the performance of the model to simulate forest regrowth after harvests is limited by the absence of regeneration and understory in the model.Most process-based models such as PnET-CN and TEM (Raich et al., 1991) have been mainly developed to simulate C balance for mature forests over the past decades (Landsberg, 2003), resulting in no provision for simulating regeneration such as shrub component and species succession in these models.Change in forest composition (e.g., evergreen and deciduous tree species and understory shrubs) along the course of succession might not be fully considered by most ecosystem models.PnET-CN does not simulate shrubs and herbs that likely dominate stands in the early successional stage after stand-replacing harvests.For this reason, the model is not able to simulate the particularly high estimated GPP and ER in the young forests where forest canopy has not yet fully recovered, if the flux partitioning algorithm used in these sites gave unbiased estimates of GPP and ER.
Understory layer is also an important component for mature forest ecosystems in terms of C fluxes (Misson et al., 2007) and stocks.Misson et al. (2007) reported that understory can contribute 11% (range, 0-39%) of GPP at 10 sites across a wide range of forest type and climate.PnET-CN slightly overestimated overstory LAI for the mature DBF sites and reasonably predicted foliar nitrogen concentration compared to satellitebased estimates (data not shown).The lack of understory layer in the model is possibly responsible for the underestimation of GPP for mature DBF sites.Species competition and cohort methods that have been employed in other models such as ED (Medvigy et al., 2009) and LPG-Guess (Smith et al., 2001)  canopy light attenuation constant coefficient is typically measured in mature forests (e.g., Ryu et al., 2008), although the coefficient is known to change with canopy cover (Brantley and Young, 2007).The use of the canopy light attenuation coefficient measured in mature forest for whole forest life simulations could underpredict stand development due to the underestimation of photosynthesis in young forests.Understanding the relationship between such parameters and state variables (e.g., LAI) is thus one of the challenges to simulate the effects of stand-replacing harvests on forest C dynamics.Third, changing climate conditions can also affect the values of some parameters.For example, wood turnover rate (%, tree mortality in terms of biomass losses) to which wood living biomass C and soil organic C are sensitive, could be altered by extreme weather conditions including droughts (Allen et al., 2010;Wang et al., 2012a).Most process-based models are not able to simulate some mechanistic processes (e.g., tree mortality and phenology) under changing climate conditions (McDowell, 2011;Richardson et al., 2011;Wang et al., 2012a), although interest in mechanically modeling forest mortality is rapidly growing (e.g., McDowell et al., 2013;Powell et al., 2013).Recent studies have revealed that climate and disturbance legacies govern forest C dynamics (Magnani et al., 2007;Bond-Lamberty et al., 2013).Future modeling efforts can benefit from improved understanding effects of climate change on the values of parameters that are assumed to be constant in the model.
Finally, the whole silvicultural system (e.g., harvests) is not fully considered in the model.Harvest methods depend on forest types, management needs, and species to be regenerated.For example, selective harvesting or shelterwood system normally is used for hardwoods in Wisconsin (Wisconsin Department of Natural Resources, 2011).Stand-replacing harvesting was assumed for both DBF and ENF chronosequence sites due to the lack of harvesting information and the types of clearing applied to the sites studied.Sensitivity analysis conducted in this study suggests that harvest intensity affects C dynamics in early succession after harvesting.Observations in residuals and post stands after each operation type (e.g., pre-commercial thinning and selective har- Full vesting) are needed to parameterize process-based models for better process-based understanding of the effects of forest harvests.

Conclusions
The PnET-CN model is generally able to simulate the effects of stand-replacing harvests on forest C dynamics (C fluxes and AGC) for two northern temperate forest chronosequences.The predicted dynamics in NEP and AGC due to clearcuts generally follow the hypothesized trajectories (Odum, 1969;Chapin et al., 2002) Full      1910 1920 1930 1940 1950 1960 1970 1980 1990 1).The time series started from the earliest major disturbance for each site.Symbols represent measured LAI.
Discussion Paper | Discussion Paper | Discussion Paper | Abstract Stand-replacing disturbances including harvests have substantial impacts on forest carbon (C) fluxes and stocks.The quantification and simulation of these effects is essential for better understanding forest C dynamics and informing forest management in the context of global change.We evaluated the process-based forest ecosystem model, PnET-CN, for how well and by what mechanisms changes of ecosystem C fluxes, aboveground C stocks (AGC), and leaf area index (LAI) arise after clearcuts.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | indicates a good agreement and a value of 100% suggests a poor agreement.The Willmott index of agreement (d ) is an indicator of modeling efficiency and is expressed Discussion Paper | Discussion Paper | Discussion Paper | . The Willmott index of agreement between simulated and tower C fluxes for both plant function types ranged from 0.91 to 0.94 with the exception of NEP (d = 0.73, n = 235).The model underestimated GPP for the DBF sites and predicted ER fairly well for all DBF sites except for the intermediate-aged hardwood site, IHW.As a result, the model underestimated NEP for most DBF sites.For IHW, the model substantially underestimated both GPP and ER but predicted NEP fairly well.For the ENF sites, the model underestimated GPP.The model predicted ER fairly well for YRP (eight years old), YJP (15-16 years old) and IRP (23 years old), but overestimated ER for the older MRP sites.Thus, the model underestimated NEP for the ENF sites.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | after fire, and then generally leveled off at 5.3-7.2m 2 m −2 from 23 to 154 years after the disturbance.A modeling study based on a modified version of Biome-BGC (Bond-Lamberty et al., 2005) also showed a similar successional change in LAI for boreal DBF and ENF.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ters et al.
could be used to improve the regeneration and understory components of PnET-CN in the future.Second, parameter values used in the model were generally derived from specific measurements for a given stand development stage particularly mature forests, although the parameter values likely differ with stand development.For example, the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | comments to the early draft of the paper.We also thank Lucie Lepine, Zaixing Zhou, Andrew Ouimette, and Alexandra Thorn for helpful discussion.Discussion Paper | Discussion Paper | Discussion Paper |

a
Estimated year of disturbance based on Ameriflux site description in AmericFlux.b Sum of wood and foliage biomass carbon from Curtis et al. (2002).c Estimated values based on measurements in 1998 to 2000 and 2002 from Cook et al. (2008).d Value in 2003 from Gough et al. (2008).e Calculated based on multi-year (1999-2003) estimations with litter traps from Gough et al. (2008Discussion Paper | Discussion Paper | Discussion Paper |

Figure 6 .
Figure 6.Simulated trajectories of LAI for each site based on the site disturbance history (Table1).The time series started from the earliest major disturbance for each site.Symbols represent measured LAI.

The Supplement related to this article is available online at doi:10.5194/bgd-11-8789-2014-supplement. Introduction
, though our simulations show that the decline in NEP is due to relatively flatting GPP and gradually increasing ER.Furthermore, our study shows that DBF recovered faster (11 years) from net C sources to net sinks and accumulated more C in AGC than ENF.Northern temperate ENF are more vulnerable to stand-replacing harvests than northern temperate DBF.Future research is needed to understand how respiration components (e.g., ecosystem and soil respiration) and production components (e.g., overstory and understory) change with forest age and their determinants, providing useful information for process-based models.Modeling the combined effects of climate change and forest management requires incorporation of forest population dynamics (e.g., regeneration and mortality), relationships between age-related model parameters and state variables (e.g., LAI), and silvicultural system into the model.Process-based ecosystem models thus could predict more reliable regional C balance under changing climate associated with disturbance regime.

Table 1 .
Site characteristics for two chronosequences of deciduous broadleaf forests (DBF) and evergreen needleleaf forests (ENF) in Upper Midwest region of Wisconsin and Michigan, USA.

Table 2 .
PnET model performance in monthly carbon fluxes (GPP: gross primary productivity; ER: ecosystem respiration; NEP: net ecosystem productivity), leaf area index (LAI), and aboveground carbon stock (AGC) for the two chronosequences.
a Normalized root mean square error.b Willmott index.c Mean value of evaluating statistics for all tested variables.Figures Back Close Full