Quantifying the effects of harvesting on carbon fluxes and stocks in northern temperate forests

Harvest disturbance has substantial impacts on forest carbon (C) fluxes and stocks. The quantification of these effects is essential for the better understanding of forest C dynamics and informing forest management in the context of global change. We used a process-based forest ecosystem model, PnET-CN, to evaluate how, and by what mechanisms, clear-cuts alter ecosystem C fluxes, aboveground C stocks (AGC), and leaf area index (LAI) in northern temperate forests. We compared C fluxes and stocks predicted by the model and observed at 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 normalized root mean square error (NRMSE) and the Willmott index of agreement (d) for carbon fluxes, LAI, and AGC in the two chronosequences were 20 % and 0.90, respectively. Simulated gross primary productivity (GPP) increased with stand age, reaching a maximum (1200–1500 g C m yr) at 11–30 years of age, and leveled off thereafter (900–1000 g C m yr). Simulated ecosystem respiration (ER) for both plant functional types (PFTs) was initially as high as 700–1000 g C m yr in the first or second year after harvesting, decreased with age (400– 800 g C m yr) before canopy closure at 10–25 years of age, and increased to 800–900 g C m yr with stand development after canopy recovery. Simulated net ecosystem productivity (NEP) for both PFTs was initially negative, with net C losses of 400–700 g C m yr for 6–17 years after clear-cuts, reaching peak values of 400–600 g C m yr at 14–29 years of age, and eventually stabilizing in mature forests (> 60 years old), with a weak C sink (100– 200 g C m yr). The decline of NEP with age was caused by the relative flattening of GPP and gradual increase of ER. ENF recovered more slowly from a net C source to a net sink, and lost more C than DBF. This suggests that in general ENF may be slower to recover to full C assimilation capacity after stand-replacing harvests, arising from the slower development of photosynthesis with stand age. Our model results indicated that increased harvesting intensity would delay the recovery of NEP after clear-cuts, but this had little effect on C dynamics during late succession. Future modeling studies of disturbance effects will benefit from the incorporation of forest population dynamics (e.g., regeneration and mortality) and relationships between age-related model parameters and state variables (e.g., LAI) into the model.


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).Harvesting is an important anthropogenic disturbance shaping North American forest Published by Copernicus Publications on behalf of the European Geosciences Union.
landscapes.Approximately 61 000 km 2 of forests were affected by harvests every year during the 2000s (Masek et al., 2011).Harvests alter forest age structure and the forest carbon (C) cycle (Magnani et al., 2007;Pan et al., 2011;Williams et al., 2012;Zhou et al., 2013a).Quantifying the effects of harvest disturbances under the context of climate change is essential for predicting forest C dynamics, informing climate policy-making, and improving forest management.Here, we focus on the assessment of an ecosystem model's ability to assess the C cycle response to clear-cut harvesting.
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 its maximum in middle age, and declines slightly thereafter in response to nutrient limitations and aging (Odum, 1969;Chapin et al., 2002;Tang et al., 2014).The successional change in plant respiration (i.e., autotrophic respiration) after stand-replacing harvesting is similar to that of GPP, although the C use efficiency, the ratio of net primary productivity (NPP) to GPP (or NPP / GPP), generally declines with forest age (DeLucia et al., 2007).As a result of these patterns, along with increases in C loss through woody litterfall, living tree biomass C gradually increases following a typical logistic growth curve (Odum, 1969;Sprugel, 1985).
Heterotrophic respiration following stand-replacing harvesting can 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 an increase in litter quantity depending on harvest types (e.g., stem-only harvesting).Heterotrophic respiration is expected to gradually decrease thereafter because the regrowing forest reduces net radiation, 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, however, heterotrophic respiration could be enhanced 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) can also be influenced by harvest types and forest composition.Unlike GPP or NPP, quantifying the trajectory of heterotrophic respiration (and consequently total ecosystem respiration) with age is not as straightforward (Amiro et al., 2010).Observational studies have shown that forest ecosystems generally become C sources (i.e., negative net ecosystem productivity, NEP) immediately following standreplacing 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 Chapin et al. (2002).
The changes in C fluxes and stocks after harvesting have been examined in many forest ecosystems using ecological measurements (e.g., eddy covariance or EC observations) from chronosequences using a space-for-time substitution approach (e.g., Gough et al., 2007;Goulden et al., 2011).The trajectory 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 3 years 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 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 into more process-based and quantitative understanding of disturbance effects on the C cycle using ecosystem models (Goulden et al., 2011).Process models require evaluation of 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 the lack of biological and climatic data (Yanai et al., 2003;Bond-Lamberty et al., 2006) and 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 clear-cuts 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 stand-replacing disturbances for forest chronosequences.
The objectives of this study were to evaluate the ability of an ecosystem model to capture the trajectories of forest C dynamics after stand-replacing harvests for two northern temperate plant functional types (PFTs: deciduous broadleaf forests, DBF; evergreen needleleaf forests, ENF), to examine which processes influence successional trajectories in these ecosystems and to test the role of PFT on the successional trajectory of C fluxes.We applied a process-based forest ecosystem model, PnET-CN (Aber et al., 1997;Ollinger et al., 2002), for simulating the effects of clear-cut on forest C dynamics, and evaluated the simulated C fluxes and stocks for both PFTs using in situ measurements (e.g., EC observations and aboveground biomass C, AGC).We hypothesized that (1) both DBF and ENF will have similar successional patterns in C fluxes (GPP, ER, and NEP) and aboveground biomass C stocks 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 field data
Our study sites consist of eight EC sites in the Upper Midwest region of northern Wisconsin and Michigan (Chen et al., 2008; Table 1).The study area is characterized by a humid-continental 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 are less than 100 years old in this region, having regenerated following past harvesting (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.Dominant tree species are bigtooth aspen (Populus grandidentata), trembling aspen (Populus tremuloides), sugar maple (Acer saccharum), red maple (Acer rubrum), red oak (Quercus rubra), basswood (Tilia american), and green ash (Fraxinus pennsylvanica).The four ENF sites also represent a chronosequence with stand age ranging from 8 to 66 years.Red pine and jack pine (Pinus banksiana) are the dominant tree species in the four ENF sites.At the two chronosequences, most sites were initiated by stand-replacing harvests.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 EC flux tower sites (Table 1).Harmonized level 4 data were used in this study.These flux data have been described and used in our previous studies (e.g., Noormets et al., 2007;Chen et al., 2008;Desai et al., 2008;Xiao et al., 2011Xiao et al., , 2014)).We also obtained leaf are index (LAI) and AGC data from the literature for each site (Table 1).

Model description
The PnET-CN model is a process-based forest ecosystem model designed to simulate C, nitrogen (N), and water dynamics at daily to monthly time steps.PnET-CN is driven by climate variables (temperature, precipitation, and photosynthetically active radiation; PAR), site variables and atmospheric properties (soil moisture, disturbance history, wet and dry N deposition, and atmospheric CO 2 concentration), and vegetation input parameters describing physiological and structural plant traits (Aber and Driscoll, 1997;Aber et al., 1997;Ollinger et al., 2002).The model has been applied and tested in the USA and Europe for simulating the 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).
A characteristic feature of PnET-CN is its use of generalized leaf trait relationships to simulate potential photosynthesis in a multilayered canopy (Aber and Federer, 1992).Actual photosynthesis is then constrained by air temperature, vapor pressure deficit, and soil water availability for simulating actual GPP.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, wood, fine roots, woody debris, and soil organic matter) and two of which are non-structural C pools stored in wood and fine roots.Photosynthetic production is allocated to each living plant component (i.e., foliage, wood, and fine roots) and to growth and maintenance respiration.Living biomass is transferred to dead woody biomass and/or to soil organic C through leaf, root and wood 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.
PnET-CN includes a complete N cycle, and simulates N mineralization and nitrification, plant N uptake, allocation, and leaching losses.N deposition is added to corresponding soil N pools (NH 4 and NO 3 ).As with C pools, N is divided into five structural pools (foliage, wood, fine roots, woody debris, and soil organic matter) and one non-structural N pool stored in the trees.C and N cycles interact closely in the model.High leaf N concentration increases net photosynthesis rate in the absence of water stress, thereby resulting in the high demand for non-structural N in plant tissues (Aber et al., 1997).When plant non-structural N is low, plant N uptake efficiency from available soil mineral N 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.
The model also simulates key hydrological processes including rainfall interception, evaporation, transpiration, surface runoff, and drainage at each time step.Rainfall  -year (1999-2003) estimations with litter traps from Gough et al. (2008).
interception is treated as a constant fraction of precipitation.
Transpiration is estimated based on water use efficiency and plant demand via photosynthesis.Surface runoff is calculated 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 (live and dead), and the loss rate of soil organic matter.In this study, when stand-replacing disturbance events occur, a uniform PFT was assumed to be regenerated on-site.For the first year after clear-cuts, a minimum LAI of 0.1 was assumed to regulate maximum potential foliage mass that controls leaf production.The photosynthetic production is transported to plant non-structural C pool where C could be allocated to leaves, stems, and roots.There is, therefore, 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).

Model inputs
The model inputs include air temperature, precipitation, PAR, wet and dry N 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 air 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 radia-tion (RAD, MJ m −2 day −1 ) using the empirical relationship (PAR = 2.05 RAD) (Aber et al., 1996).The data from 1981 through 2010 were repeated as needed to create the time series from 1850 to 1980.
Annual rates of wet and dry N deposition were obtained from the United States Environmental Protection Agency (EPA; http://java.epa.gov/castnet/clearsession.do).The N deposition rates were measured at the Wellston station (44.22 • N, 85.82 • W) for the period 1994-2011.We also obtained the N 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 N deposition rates.Monthly wet deposition rates, needed for the model, were generated from annual wet N deposition through the weighted coefficients (the ratio of monthly precipitation to total precipitation from March to November).We assumed that there is no wet N deposition in the winter.The soil water holding capacity in the rooting zone (100 cm) for each site was derived from the gridded multilayer soil characteristics data set (STATSGO, Miller and White, 1998).For the period 1959-2010, we used the CO 2 concentrations data from Mauna Loa.For the time period 1901-1958, we derived the time series of the historical atmospheric CO 2 mixing ratio using a spline fit to the ice-core record (Etheridge et al., 1996), as described by McGuire et al. (2001) and used by Xiao et al. (2009).We used the CO 2 concentration in 1901 for the simulation period prior to 1901.
For each site, we prescribed the disturbance events using the site disturbance history (Table 1).For each standreplacing 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.The soil removal fraction was assumed to be 0, given that the content of soil organic C might not be considerably affected by harvesting (Johnson and Curtis, 2001;Yanai et al., 2003).We also conducted a sensitivity analysis on these assumptions as described below in Sect.2.4.

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 N mineralization, and foliar N concentrations.The parameter values used in this study are given in Table S1 in the Supplement.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.Thirty-year climate normals (1981-2010), preindustry N deposition rates, and historical CO 2 concentrations (1901) were used for the spin-up runs.
To examine the effects of stand-replacing harvests on C dynamics, we conducted all simulations using the site disturbance history (Table 1), vegetation parameters (Table S1), climate, N 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 EC flux sites.We used two statistical measures to evaluate the overall model performance: the normalized root mean square error (NRMSE) and the Willmott in-  dex of agreement.The NRMSE (Eq. 1) was used to assess the difference between predicted (P ) and observed (O) variables, and can be expressed as 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 indicates perfect agreement and a value of 100 % suggests poor agreement.The Willmott index of agreement (d) is an indicator of modeling efficiency and is expressed as (2) A value of 1 indicates perfect agreement, and a value near 0 indicates weak agreement (Willmott, 1982).The sensitivity of ecosystem C dynamics to changes in harvesting practices during the secondary succession was assessed using sensitivity analysis.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) and soil organic matter loss (20 and 40 %, compare to 0 in the model test) to assess the effects of different harvest intensity and soil organic matter loss on C dynamics.We also tested the model sensitivity to CO 2 fertilization for evaluating potential climate change effects.

Evaluation of modeled carbon fluxes and stocks
Simulated C fluxes were generally consistent with EC derived C fluxes for both DBF and ENF sites (Figs. 1 and 2).The NRMSEs between simulated and tower fluxes (GPP, ER, and NEP) were between 10 and 21 % (Table 2).The Willmott index of agreement between simulated and tower C fluxes for both PFTs 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 (8 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.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.The NRMSE was 28 for AGC and 31 % for LAI.The Willmott index of agreement was 0.95 and 0.96 for AGC and LAI, respectively.Overall, the model evaluation metrics indicated that the model performed better in the DBF sites than in the ENF sites.

Effects of clear-cutting on carbon fluxes and stocks
PnET-CN generally captured the changes of C fluxes following the clear-cuts for each chronosequence site (Fig. 4).
The predicted annual GPP generally increased with time since disturbance and approached 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 age older than 60 years, the predicted annual ER for both PFTs showed a relatively  1).The time series started from the earliest major disturbance for each site.
flat pattern, contrary to theoretical expectations, arising from the small amount of change of both autotrophic and heterotrophic respiration with age (Fig. S1 in the Supplement).
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 clear-cutting for both DBF and ENF (Fig. 5).Within ∼ 6 years for the DBF sites and ∼ 17 years for the ENF sites, the GPP : ER ratio gradually increased and its average value became larger than 1 (NEP > 0).The simulated peak GPP : ER ratio for DBF (1.6) occurred 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 first 6 and 17 years after stand-replacing 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.
LAI fully recovered within 10-15 years after disturbance for the DBF sites and within 40 years of age for the ENF sites (Fig. 6).The 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 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 analysis
Harvest intensity had little effect on long-term C dynamics for both PFTs, but it had sizeable effects during early succession (Fig. 8).Increasing harvest intensity delayed GPP (Fig. 8a and f) and LAI (Fig. 8d and i) rises and led to a lower reduction in ER (Fig. 8b and g), resulting in laterrising NEP (Fig. 8c and g 2010) using eddy covariance observations from forest chronosequences in North America.Solid and hollow circles represent measured annual and growing season (May to October) ratios, respectively.The simulated curves were smoothed using a moving average filter with a span of 5.
100 % removal of living trees) also directly reduced living tree AGC (Fig. 8e and i).By reducing the 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 an increase in AGC about 12 and 16 %, respectively, after a 100-year harvest cycle.For ENF, average annual NEP decreased about 1 %, and AGC decreased nearly 6 % for both reduced harvest intensity scenarios.Increasing the soil removal fraction parameter resulted in lower GPP and ER along succession and lower NEP in middle succession for both DBF and ENF (Fig. S2a-c and f-h in the Supplement).Greater soil removal fraction enhanced the leaf biomass of DBF in middle and late succession (Fig. S2d), but restricted the leaf biomass of ENF in late succession (Fig. S2i).Increasing the soil removal fraction parameter (20 and 40 % removal of soil organic matter) strongly reduced living AGC (16 and 39 %, respectively) for DBF (Fig. S2e) but slightly decreased living AGC (up to 5 %) for ENF (Fig. S2j).There were insignificant effects of CO 2 fertilization on carbon dynamics for both DBF and ENF in our sensitivity analysis (Fig. S3 in the Supplement).

Carbon fluxes and stocks following clear-cutting
PnET-CN generally simulated the expected post-harvest trajectories in C fluxes (GPP, ER, and NEP) and stock (LAI and AGC).Our simulations showed that LAI first increased rapidly and then stabilized during the following development stages, because the model estimates foliage growth through the parameter of maximum relative growth rate (Table S1) with the restriction of current foliage biomass and resource availability.This modeled response is consistent with the previous finding that foliage biomass increased rapidly after disturbance and then stabilized (Sprugel, 1985).Our chronosequence-based results are generally consistent with previous results.For example, Goulden et al. (2011) observed that LAI along a chronosequence of boreal forest stands increased rapidly from 0.3 m 2 m −2 1 year after fire, and then generally leveled off at 5. The simulated successional change in annual GPP for both PFTs generally followed the trajectory hypothesized by Odum (1969).However, despite a slight decrease in GPP hypothesized in Odum's trajectory, 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 also stabilizes.Our results are consistent with previous studies showing relatively flat patterns in GPP 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 jack pine forests in Saskatchewan (Zha et al., 2009).
Furthermore, Humphreys et al. (2006) reported continuous increases of GPP with increasing forest age for temperate rainforests using three different stands at different stages of development (2, 14, and 53 years of age) following clear-cuts in British Columbia, Canada.However, northern temperate ENF showed a small difference in GPP between young and mature sites (Noormets et al., 2007;Desai et al., 2008).Desai et al. (2005) found that a nearby old-growth mixed forest had slightly lower GPP and significantly higher ER than nearby DBF sites.Site-to-site variations in species and soil fertility could result in variations in the successional trajectory of GPP after clear-cuts such that the observed trajectories  1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 1).The time series started from the earliest major disturbance for each site.Symbols represent measured LAI.
may deviate from hypothesized or modeled trajectories.But the model was unable to simulate high GPP rates estimated by the EC technique in mature forests regardless of vegetation type, suggesting that there is room for improvement in model simulation of secondary succession.In addition, our chronosequences lack old-growth sites and do not encompass the full range of forest development stages, which limits the representativeness of the C flux and stock trajectories derived from chronosequence studies based on EC or other ecological observations (e.g., Clark et al., 2004;Humphreys et al., 2006;Noormets et al., 2007).We found that annual ER for secondary temperate forests declined slightly in the first ten years because of low autotrophic respiration at first after the removal of trees.Amiro et al. (2010) reported that ER was reduced in the very first year following harvests for a number of EC flux sites over North America.Previous field studies showed that ER following clear-cuts increased with forest age (e.g., Humphreys et al., 2006;Zha et al., 2009), partly supporting our results that ER slightly increased after the short decline period (10-25 years of age) in northern temperate forests until the stands reached maturity.Martin and Bolstad (2005) showed 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 to 2006 in the same region.Tang et al. (2009) reported that growing-season soil respiration in 2005 was 690 g C m −2 in a mature DBF near WIC tower site, which is in the range of our simulations (Fig. S1).The model underestimated GPP, but estimated ER well for mature DBF sites, indicating that the model likely overestimated autotrophic respiration.EC derived ER were  1910 1920 1930 1940 1950 1960 1970 1980 1990 2000 1).
The time series started from the earliest major disturbance for each site.Symbols represent estimated AGC.
usually lower than chamber-based estimates at 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 overestimated ER for the mature site because of an overestimated soil decomposition rate.Our simulations also show that DBF had a slightly higher soil respiration rate than ENF (Fig. S1), which is consistent with the finding that chamber-based soil respiration was slightly higher for DBF than for ENF in Wisconsin (Euskirchen et al., 2003).The changes of ER in secondary forests after clear-cutting differ among sites because of different site conditions (e.g., quantity and quality of soil organic C and litter C) 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 fire and harvest 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 underestimated the observed ratios for mature sites.Growing season GPP : ER ratios are typically 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 for the EC flux data may lead to the overestimation of net ecosystem exchange due to the lack of winter C flux observations for the ENF sites and two of the DBF sites (YHW and IHW).(f-j) evergreen needleleaf forests (ENF) at a mature red pine site over a 100-year harvest cycle.The simulated curves were smoothed using a moving average filter with a span of 5.
Our simulated successional dynamics of NEP following clear-cuts generally supported the trajectories of Chapin et al. (2002), but for different reasons.The hypothesized trajectories show declining GPP and relatively flat ER with time.Our simulated decline in NEP resulted 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).A recent North American Carbon Program (NACP) synthesis study showed similar changes in NEP after either stand-replacing fire or harvest based on EC chronosequence measurements across North America (Amiro et al., 2010).Chapin et al. (2002) hypothesized that heterotrophic respiration is initially high following disturbance, declines in early succession, rises in middle succession, and declines thereafter, while NPP reaches a peak in middle age and declines in old stands.The simulated successional trajectories in heterotrophic respiration were supported by hypothesized pattern change, although our simulated NPP did not decline in mature stands (Fig. S1).Previous studies also 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.A decline of NPP with age was not predicted in this study.Recent studies found that the decline of NPP with age is not "universal" (Ollinger and Smith, 2005;He et al., 2012).Validation of the simulated NPP was not possible in this study due to the lack of NPP measurements 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 that NEP continually increased with age because of 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 similar regions.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 that 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 sites, including a northern hardwood forest of central Massachusetts (200 ± 40 g C m −2 yr −1 , Barford et al., 2001) and four eastern North American DBF sites (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 showed 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.However, the representativeness and generalization of these findings were limited by the small number of young stands (Powers et al., 2012).
Sensitivity analysis shows that more intensive harvests could have larger and longer impacts on successional trajectories of C dynamics in early succession for both DBF and ENF.Fewer EC-based studies have investigated the effects of harvest intensity on forest C fluxes (e.g., GPP, ER, and NEP) because of the high establishment cost of EC systems.Nevertheless, some modeling studies have provided insights into how forest C fluxes and stocks are affected by harvest intensity.Our findings are supported by previous modeling studies.For example, a recent modeling study of temperate forests reported that more intensive harvests increased the recovery time of NPP for ENF and DBF in Minnesota and Wisconsin, USA (Peters et al., 2013).In the boreal forests of central Canada, less intensive harvest and longer rotation length might increase total C sink up to 40 % (Peng et al., 2002), although recent studies indicate that longer rotation length could not necessarily increase C sequestration under changing climate conditions (Wang et al., 2012b;Wang et al., 2013).A recent synthesis study on the effects of partial cutting on forest carbon stocks found that partial cutting has no significant effects on litter C and soil organic C, although more intensive cutting can reduce AGC more (Zhou et al., 2013b).This synthesis study is not able to determine the recovery duration due to the lack of long-term observations.If harvesting operations largely reduce soil organic matter, C fluxes (e.g., GPP, NPP, ER, and NEP) and living AGC are reduced for both PFTs.Consistent with this, Peters et al. (2013) showed that simulated NPP could not recover to pre-harvesting levels due to greater removal of soil organic matter.Therefore, our model results suggest management practices should aim to decrease soil disturbance caused by harvest operations.

Differences between deciduous broadleaf and evergreen needleleaf forests
We found that DBF may reach a peak in LAI and GPP faster than ENF after clear-cutting, showing clear differences in pattern of ecosystem development between the DBF and ENF sites.More rapid recovery of LAI and GPP for DBF sites lead to sooner recovery of NEP and AGC regardless of harvest intensity, supporting our second hypothesis.The foliage-related parameters such as FolRelGroMax and AmaxB mainly govern the differences in successional trajectories between the two PFTs (Table S1).DBF is assumed to have more productive foliage than ENF, and more photosynthetic production then can 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 studies showing that the temperate DBF in northern Michigan rapidly became a net C sink after 6 years following disturbances (Gough et al., 2007) and that ENF stands in northern Wisconsin became net C sinks within 10-15 years after harvesting (Noormets et al., 2007).Through the analysis of the Forest Inventory and Analysis (FIA) data, Williams et al. (2012) suggested that faster growth in AGC at high pro-ductivity sites caused higher C fluxes and stocks.Our findings are also consistent with a recent modeling study suggesting that temperate DBF switches to positive NEP faster than temperate ENF after clear-cuts, 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 more years to attain a positive NEP than boreal DBF after clear-cuts in Saskatchewan, Canada (Wang et al., 2012b).These observed and modeled successional changes further indicate that DBF tend to have higher photosynthetic capacity than ENF in the early stage of stand development following stand-replacing harvests.
The sensitivity analysis suggests that more productive forests could be more strongly affected by greater soil removal fractions, as soil removal reduces soil organic matter thereby resulting in relatively low N mineralization in the model.Peters et al (2013) showed that NPP was more strongly reduced for aspen than for jack pine in their simulations.However, productive lodgepole pine (Pinus contorta Dougl.ex Loud.) could maintain high productivity at 12 years after harvest disturbance regardless of soil organic carbon removal and soil compaction treatments because beneficial ectomycorrhizal fungi associated with lodgepole pine could help access N from organic matter; hybrid white spruce (Picea glauca × engelmannii [Moench] Voss), meanwhile, was more sensitive to the treatments (Kranabetter et al., 2006).The discrepancy might be caused by the lack of representation of a relationship between fungi and plants in the model.

Limitations and challenges
PnET-CN can explicitly simulate the effects of disturbance, pollution, and climate change on forest C dynamics (e.g., Ollinger et al., 2002;Pan et al., 2009;Peters et al., 2013).Despite the capability of the model, 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.
The performance of the model to simulate forest regrowth after harvests is limited by the absence of population and community dynamics associated with regeneration.Most process-based models such as PnET-CN and TEM (Raich et al., 1991) have been developed primarily to simulate C balance for mature forests over past decades (Landsberg, 2003), resulting in no provision for simulating regeneration such as shrub component and species-specific successional dynamics in these models.Changes in forest composition (e.g., evergreen and deciduous tree species and understory shrubs) along the course of succession are not 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.Understory layer is also an important component for mature forest ecosystems in terms of C fluxes 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 types and climates.PnET-CN slightly overestimated overstory LAI for the mature DBF sites and reasonably predicted foliar N concentration compared to satellite-based estimates.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) could be used to improve the regeneration and understory components of PnET-CN in the future.
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 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 slow down 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.
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 the mechanistic processes associated with tree mortality under changing climate conditions (Mc-Dowell, 2011;Wang et al., 2012a), although there is growing interest in the mechanistic modeling of forest mortality (e.g., McDowell et al., 2013;Powell et al., 2013).Recent studies have revealed that climate and disturbances govern forest C dynamics (Magnani et al., 2007;Bond-Lamberty et al., 2013).Future modeling efforts can benefit from improved understanding of the effects of climate change on parameter values that are assumed to be constant in the model.
Harvest methods depend on forest types, management needs, and species to be regenerated.For example, selective harvesting or the shelterwood system is typically 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.The 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., precommercial thinning and selective harvesting) are needed to parameterize process-based models for better mechanistic understanding of the harvest effects on forest C dynamics.

Conclusions
The PnET-CN model was 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 following clear-cuts generally follow the hypothesized trajectories, although our simulations show that the decline in NEP was due to relatively stable GPP and gradually increasing ER.Our study also 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 is more vulnerable to stand-replacing harvests than northern temperate DBF.Future research is needed to better 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.Modeling the combined effects of climate change and forest management will benefit from the 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.With these improvements, process-based ecosystem models can better simulate regional C balance associated with disturbance regime under changing climate.
The Supplement related to this article is available online at doi:10.5194/bg-11-6667-2014-supplement.

Figure 3 .
Figure 3. Comparisons of simulated and observed (a) leaf area index (LAI) and (b) aboveground carbon stock (AGC) for all eight sites.

Figure 4 .
Figure 4. Simulated trajectories of GPP, ER, and NEP for each site based on the site disturbance history (Table1).The time series started from the earliest major disturbance for each site.

Figure 5 .
Figure 5. Simulated trajectories of the annual GPP / ER ratio with stand age for (a) deciduous broadleaf forests (DBF) and (b) evergreen needleleaf forests (ENF).The dashed line is a fitted curve derived by Amiro et al. (2010) using eddy covariance observations from forest chronosequences in North America.Solid and hollow circles represent measured annual and growing season (May to October) ratios, respectively.The simulated curves were smoothed using a moving average filter with a span of 5.
3-7.2 m 2 m −2 from 23 to 154 years after stand-replacing crown fire.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.

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.

Figure 7 .
Figure 7. Simulated trajectories of aboveground biomass carbon (AGC) for each site based on the site disturbance history (Table1).The time series started from the earliest major disturbance for each site.Symbols represent estimated AGC.

Figure 8 .
Figure8.Sensitivity of carbon fluxes (GPP, gross primary production; ER, ecosystem respiration; NEP, net ecosystem production) and stand characteristics (LAI: leaf area index; AGC: aboveground carbon stock) to changes in harvest intensity (reduced by 0.2 and 0.4 compared to 1 for assumed clear-cuts used in the model tests) for (a-e) deciduous broadleaf forests (DBF) at Willow creek and (f-j) evergreen needleleaf forests (ENF) at a mature red pine site over a 100-year harvest cycle.The simulated curves were smoothed using a moving average filter with a span of 5.

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

Table 2 .
PnET-CN 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.

www.biogeosciences.net/11/6667/2014/ Biogeosciences, 11, 6667-6682, 2014 the
The model thus is not able to simulate the particularly high GPP and ER in young forests where forest canopy has not yet fully recovered.