Net ecosystem productivity of temperate and boreal forests after clearcutting—a Fluxnet-Canada measurement and modelling synthesis

Clearcutting strongly affects subsequent forest net ecosystem productivity (NEP). Hypotheses for ecological controls on NEP in the ecosystem model ecosys were tested with CO2 fluxes measured by eddy covariance (EC) in three postclearcut conifer chronosequences in different ecological zones across Canada. In the model, microbial colonization of postharvest fine and woody debris drove heterotrophic respiration (Rh), and hence decomposition, microbial growth, N mineralization and asymbiotic N2 fixation. These processes controlled root N uptake, and thereby CO2 fixation in regrowing vegetation. Interactions among soil and plant processes allowed the model to simulate hourly CO2 fluxes and annual NEP within the uncertainty of EC measurements from 2003 to 2007 over forest stands from 1 to 80 yr of age in all three chronosequences without site- or species-specific parameterization. The model was then used to study the impacts of increasing harvest removals on subsequent C stocks at one of the chronosequence sites. Model results indicated that increasing harvest removals would hasten recovery of NEP during the first 30 yr after clearcutting, but would reduce ecosystem C stocks by about 15% of the increased removals at the end of an 80-yr harvest cycle.


Introduction
Net biome productivity (NBP) of forest ecosystems can be adversely affected by an increased frequency and/or intensity of disturbances such as harvesting, fire and insects. The intensity of harvesting may vary from minimal, as in selective logging where the forest stand remains mostly in place, to severe, as in clearcut logging where the forest stand is completely replaced. Clearcut logging adversely affects NBP through: (1) direct export of C as harvested wood products (2) increased heterotrophic respiration (R h ) driven by decomposition of detritus, including above-and below-ground fine litter and woody debris (WD), left in large amounts after clearcutting and warmed by exposure to solar radiation and (3) slow recovery of net primary productivity (NPP) caused by small leaf areas and root lengths of plants regrow-ing from seeds or rootstocks under limited nutrient availability. Rates of R h may exceed those of NPP for several years after clearcutting, during which time substantial losses of ecosystem C stocks may occur. These losses may continue for 10 yr or longer in Canadian forests (Kurz and Apps, 1999;Litvak et al., 2003), for at least 14 yr in Siberian pine forests (Schulze et al., 1999) and for 14 (Janisch and Harmon, 2002) or 20 (Cohen et al., 1996) years in coniferous forests of the Pacific Northwest. Eventually R h declines as detritus remaining from the previous forest is depleted, while NPP rises as forest leaf area and root length recover, so that the forest starts to regain C.
The time course of C loss and recovery after disturbance strongly affects NBP over the forest life cycle, so that much effort has been invested in measuring and modelling rates of litter decomposition and tree growth after stand-replacing harvests or fires. Approaches to estimating WD losses have differed in complexity. Most researchers (e.g. Janisch and Harmon, 2002;Shorohova et al., 2008;Melin et al., 2009) have derived firstorder decomposition rate constants to estimate C losses from WD stocks following harvest. Shorohova et al. (2008) further proposed different rate constants for bark and wood decomposition while accounting for wood fragmentation. However other researchers (e.g. Mäkinen et al., 2006;Montes and Cañellas, 2006) have found that WD decomposition rates could better be fitted with Gomperz or Chapman-Richards functions in which a time lag was modelled before the onset of substrate-limited decomposition. These functions were based on the observation that decay is slow until heterotrophic decomposers have become established within the substrate (Grier, 1978;Harmon et al., 1986). Montes and Cañellas (2006) further parameterized these functions to allow more rapid decomposition rates in more advanced WD decay classes. Diochon et al. (2009) found that significant post-harvest C losses also occurred in the mineral soil below the forest floor, and so must be included in a full accounting of post-harvest changes in C stocks, although Johnson and Curtis (2001) found no consistent effect of harvesting on soil non-litter C.
Plant regrowth after a stand-replacing disturbance is strongly affected by mineralization of nutrients from forest litter. During the first few years after disturbance, the rapid decomposition of fine plant litter, such as foliar litter with comparatively low C:nutrient ratios can cause a transient net mineralization of nutrients that may stimulate early regrowth of LAI and rise in gross primary productivity (GPP; Kimmins, 2004), mostly by herbaceous pioneer species. However this post-harvest mineralization is sometimes not found, likely when fine litter mineralization is constrained by low pH or nutrient content (Grenon et al., 2009). The concurrent decomposition of large amounts of WD with high C:nutrient ratios immobilizes nutrients (Grier, 1978;Schimel and Firestone, 1989), limiting nutrient uptake by regrowing plants with small root systems, thereby lowering foliar N concentrations  and hence GPP (Kimmins, 2004). This constraint on nutrient uptake is gradually relieved as continued decomposition of WD generates products with smaller C:nutrient ratios, allowing eventual remineralization of WD nutrients. This constraint is relieved even further as plant roots proliferate during regrowth, improving their nutrient uptake capacity and allowing LAI and GPP to rise, mostly from dominant tree species. Thus Taylor et al. (2005) developed different Chapman-Richards functions to model the regrowth of trees, seedlings and shrubs after harvest. These functions were parameterized to allow rapid early growth and subsequent decline of seedlings and shrubs, but slow early growth and subsequent faster growth of trees. Similar functions have also been used to simulate early tree growth in inventory-driven models of forest growth (e.g. Kurz et al., 2009).
Most post-harvest models of litter decomposition and forest regrowth have been derived empirically from studies of C inventories in post-harvest chronosequences (summarized in Melin et al., 2009) so that the validity of their parameters is likely to be limited to the conditions of the studies. A more process-based model of decomposition and regrowth might enable a more ro-bust simulation of post-harvest changes in forest C stocks that could be used in life-cycle forest C studies without the need for site-specific parameterization. Such a model should be able to simulate the processes and kinetics of colonization and subsequent decomposition of new detritus stocks, as well as those of older soil organic carbon (SOC) stocks, by heterotrophic decomposers. The model should also be able to simulate the mineralization versus immobilization of N driven by this decomposition, and its effects on plant N uptake during forest regrowth.
Chronosequence studies of post-disturbance changes in C stocks can provide only weak tests of process models because these changes represent temporally aggregated results of disturbance effects on ecosystem processes, rather than effects on the processes themselves. Process-level disturbance effects may be better represented by changes in CO 2 effluxes versus CO 2 influxes measured at eddy covariance (EC) flux towers in post-disturbance chronosequences. These fluxes directly indicate rates of ecosystem respiration versus primary productivity as affected by rates of detritus decomposition versus plant regrowth, and so provide better constrained tests of process models. In this study, we test hypotheses for decomposition and regrowth in the ecosystem process model ecosys against EC CO 2 fluxes in three diverse post-clearcut chronosequences in British Columbia (BC), Saskatchewan (SK) and Quebec (QC) as part of the Canadian Carbon Program (CCP). To evaluate model sensitivity to changes in harvest intensity, the model was then used to project long-term effects of different harvest removals on forest C stocks at one site in the BC chronosequence.

General
The equations and parameters used to model decomposition and growth remain the same as those documented in several earlier modelling studies (e.g. Grant et al., 1993aGrant et al., ,b, 2007aGrant et al., ,b, 2008Grant et al., , 2009aGrant, 1998Grant, , 2004. Those equations particularly relevant to this study are described in more detail below, where key model parameters appear in bold with values given in the Appendix. A fuller description of all algorithms by which decomposition and growth are governed in ecosys may be found in 'Appendix A. Supplementary data' available online in Grant et al. (2009a).

Detritus decomposition
Decomposition of C (D SC ) from organic matter (S C ) colonized by microbial C (M C ) occurs concurrently within each of five organic matter-microbe complexes [i = woody debris, fine nonwoody litter, animal manure, particulate organic matter (POC) and humus], each with kinetic components (e.g. j = protein, carbohydrate, cellulose and lignin, where i = woody debris or litter) in each organic substrate (k = original complex matter, Tellus 62B (2010), 5 sorbed organic matter or microbial residues) in each soil or litter layer l. Decomposition (D SC ) is a temperature-dependent function (f tg ) of a specific decomposition rate (D SC ) and the active component (a) in the M C of all heterotrophic functional types (n = obligately aerobic bacteria and fungi, facultatively anaerobic bacteria, anaerobic acetogens, acetotrophic methanogens, aerobic and anaerobic diazotrophs) (first line of eq. 1). D SC is constrained by both substrate and microbial concentrations represented by half-saturation and inhibition constants (K mD and K iD , respectively) (third line of eq. 1): The latter constant inhibits D S as aqueous concentrations of M C ([M C ] = M C /θ ) rise with declining water content θ during drying of S C . Decomposition of C drives that of N and P (D SN and D SP ):

Heterotrophic respiration
Products of D SC , D SN and D SP are added to dissolved pools (Q) of which Q C is the substrate for temperature-dependent heterotrophic respiration (R h ) according to a specific respiration rate (R h ) constrained by a half-saturation constant (K mQ ) (first line of eq. 3) where R h is further constrained by microbial N and P contents (M N and M P ) with respect to their maximum set concentrations (C N and C P ), and by uptake versus demand for O 2 (U O 2 versus U O 2 ) (third line of eq. 3) as described in earlier papers cited above. This R h sustains maintenance and growth respiration (R m , a temperature-dependent function of M N , and R g , respectively) and R g drives assimilation of Q C (U C ) for growth of M C according to the energy yield ( G) versus construction energy cost for microbial growth (E n ), effectively the growth efficiency U Ci,n,l = min R hi,n,l , j R mi,n,j ,l + R gi,n,l (1 + G/E m ) (5) with associated assimilation of Q N and Q P U Ni,n,l = U Ci,n,l Q Ni,l /Q Ci,l (6a) Microbial C dynamics are given by the difference between assimilation and respiration R h (microbial growth) less first-order decay (D M ) allocated to each kinetic component j (F j ): δM Ci,n,j ,l /δt δM Ci,n,j ,l /δt = F j (U Ci,n,l − R mi,n,l ) − D MCi,n,j ,l R hi,n,l < j R mi,n,j ,l . (7b)

Mineralization and immobilization
The exchange of inorganic N or P between microbial biomass and the soil solution is determined by the ratios between M C and M N or M P resulting from U C versus U N or U P with respect to their maximum set concentrations, for example, for NH 4 and similarly for NO 3 − and PO 4 2− . Immobilization is constrained by maximum uptake rate (U NH 4 ), microbial surface area (A), minimum concentration ([NH + 4mn ]) and half-saturation constant (K NH 4 ). Thus microbial N and P dynamics are given by the difference between uptake U of organic plus inorganic N and P and first-order decay D M , that accounts for partial recycling of microbial N and P when these are limiting to microbial growth δM Ni,n,j ,l /δt = F j U Ni,nl, + U NH 4i,n,j,l + U NO 3i,n,j,l − D MNi,n,j ,l (9a) δM Pi,n,j ,l /δt = F j U Pi,nl, + U PO 4i,n,j,l − D MPi,n,j ,l . (9b)

Colonization of detritus and soil organic matter
S C added to the organic substrate (k = original complex matter) in each complex (i = woody debris, fine litter, manure, POC or humus) by litterfall or microbial transformation has to be colonized by M C before decomposition can begin. Other substrates (k = sorbed organic matter, microbial residues) in each complex i are products of microbial transformations and so are considered to be fully colonized. The colonization rate is driven by new microbial growth U C − R h in each i from eq. (7a) according to a specific colonization rate β: δS Ci,j ,k,l /δt = β n (U Ci,n,l − R hi,n,l )(S Ci,j ,k,l /S Ci,k,l ) where S C represents the remaining non-colonized k and S Ci,k,l = j S Ci,j ,k,l . The first S C term in brackets allocates colonization of S C to each kinetic component j, and the second term constrains colonization as S C declines with respect to S C according to an inhibition constant (K iS ). This rate then determines S Ci,j,k,l for k = original complex matter in eq. (1).
Colonization of new S C (e.g. litterfall) is started by microbial growth from existing S C in each complex i. S C grows as colonization proceeds, raising D SC (eq. 1) and hence R h (eq. 3), R g (eq. 4), U C (eq. 5) and U C − R h , subject to constraints on D SC by litter water content and temperature (eq. 1), and on R h by microbial temperature and nutrient uptake (eq. 3). These rises drive further growth in S C (eq. 10), causing a positive feedback cycle in which colonization, decomposition, respiration and microbial growth accelerate. However S C ' declines as colonization proceeds, gradually slowing and eventually reversing this acceleration as opportunity for further colonization diminishes (eq. 10). This algorithm thereby reproduces the dynamics of woody decomposition represented empirically by Gomperz or Chapman-Richards functions in other studies (Mäkinen et al., 2006;Montes and Cañellas, 2006). Values of β and K iS in eq. (10)) were selected to allow large deposits of WD to become fully colonized, corresponding to decay class V as defined by Siitonen et al. (2000), within 10-30 yr of deposition depending on climate, as estimated in WD decomposition models fitted to field measurements by Montes and Cañellas (2006), Müller-Using and Bartsch (2009) and Ranius et al. (2003). However these same values of β and K iS allow more rapid colonization of fine litter in which D S (eq. 1), R h (eq. 3), U C (eq. 5) and hence microbial growth U C − R h is faster.

Gross primary productivity and autotrophic respiration
Canopy GPP (CO 2 fixation) is the sum of that for each leaf surface defined by population, branch, node, layer, azimuth, inclination and exposure (sunlit or shaded), calculated from Farquhar et al. (1980) and coupled to plant water and nutrient status as described in eqs (C1)-(C11) from Grant et al. (2007b) and in eqs (A7f,g,h) in Grant et al. (2010). Products of GPP are added to nonstructural C pools σ C in branches from where they are transferred to σ C in roots and from there to mycorrhizae according to concentration gradients of σ C generated by production versus autotrophic respiration R a . Below-ground R a is the sum of that from σ C (R c ) and that from remobilization of protein structural C (R s ) by population (i = species or cohort) and organ (j = roots or mycorrhizae) in each rooted soil layer l of each order (m = primary or secondary) and node (n) where R c is a temperature-dependent function (f ta ) of σ C constrained by uptake (U) versus demand (U ) for O 2 and by comparative conductance to nonstructural C transfer g (Grant, 1998) and R s is driven from maintenance respiration requirements (R m ) unmet by R c R si,j ,l,m,n = − min{0, R ci,j ,l,m,n − R mi,j ,l,m,n }, where R s drives litterfall l C of non-protein structural C associated with protein C remobilized by R s . Growth respiration R g is the excess of R c over R m constrained by root or mycorrhizal turgor (ψ t ) above a threshold value (ψ t ) R gi,j ,l,m,n = max{0, min{(R ci,j ,l,m,n − R mi,j ,l,m,n ) × min{1, max{0, ψ ti,j ,l − ψ t }}.
2.7. Plant growth R g drives root and mycorrhizal growth according to their respective growth yields (Y g ) estimated from biochemical composition (e.g. Waring and Running, 1998) δM Ci,j ,l,m,n /δt and associated assimilation of N and P according to set ratios for roots. Growth rate (δM C /δt) drives extension of root and mycorrhizal lengths (L) according to their specific volumes (ν), densities (ρ), internal porosities (θ P ), radii (r) and an assumed cylindrical geometry δL i,j .l,m,n /δt = δM Ci,j ,l,m,n Root and mycorrhizal lengths and resulting surface areas (A) determine uptake (U) of inorganic N and P, for example, NH 4 + in eq. (17), by determining the nutrient concentrations at root and mycorrhizal surfaces ([NH + 4i,j ,l ]) at which radial transport by mass flow and diffusion, driven by the nutrient concentrations in the soil solution ([NH + 4l ]) (eq. 17a), equals active uptake by the root and mycorrhizal surfaces, driven by maximum specific uptake (U ) and the half-saturation constant (K NH 4 ) (eq. 17b) Parameters for root and mycorrhizal uptake (U NH 4i,j,l in eq. 17b) are the same as those for microbial uptake (U NH 4i,n,j,l in eq. 8b) with which U NH 4i,j,l in eq. (17b) competes. Products of U are added to nonstructural N and P pools (σ N and σ P ) in root and mycorrhizae which are coupled with σ C generated from CO 2 fixation, in mycorrhizae, roots and branches. Transfers among these pools (eq. A7a,b in Grant et al., 2010) are driven by concentration gradients generated by acquisition versus consumption of nonstructural N and C in mycorrhizae, roots and branches (eq. A7c,d in Grant et al., 2010). Ratios of nonstructural N and C in branches govern CO 2 fixation (eqs A7f,g,h in Grant et al., 2010) by (1) setting ratios of structural N and C in leaves (eq. A7e in Grant et al., 2010) and hence maximum carboxylation rates (eqs A7i,j in Grant et al., 2010), and (2) determining rubisco activation through product inhibition (eq. A7k in Grant et al., 2010). R a of above-ground phytomass is calculated from non-structural C pools for each species (i), branch (j), organ (k = leaves, twigs, branches, boles, reproductive) and node (n) in the same manner as is the R a of roots and mycorrhizae in eqs.
(11)-(14) above. Above-ground growth is thus modelled as with associated assimilation of N and P, calculated for organ k = leaf from ratios of nonstructural C, N and P as δM Pi,j ,k,n /δt from which ratios of nonstructural N and C are derived as M N :M C . Growth in leaf mass drives expansion of leaf area (A) constrained by leaf turgor (ψ t ) in eq. (20), assuming uniform growth of individual leaves (M C divided by population y) in three dimensions (Grant and Hesketh, 1992) Similarly, growth in twig, branch and bole masses drive extension of twig, branch and bole lengths (Grant and Hesketh, 1992). Values of A i,j,k,n, further resolved into layer, azimuth and inclination, are used to calculate radiation absorption and hence CO 2 fixation coupled to leaf water and nutrient status (Grant, 2004;Grant et al., 2007aGrant et al., ,b, 2009aGrant et al., ,b, 2010. Each plant population is initialized only with a small nonstructural C reserve at planting which is transferred to σ C in branches and roots during germination. These σ C then drive initial R a (eq. 12) and δM C /δt (eq. 15) and hence root elongation (eq. 16) and leaf expansion (eq. 20) required for nutrient uptake (eq. 17) and CO 2 fixation. Upon exhaustion of this reserve, each population must sustain further nutrient uptake and CO 2 fixation from root elongation and leaf expansion driven by the non-structural products of nutrient uptake and CO 2 fixation. No areas or lengths of roots, mycorrhizae or leaves are prescribed. This creates feedback during plant growth, which can be positive when growth exceeds litterfall, or negative when growth does not.

Model initialization and spinup
The different sites at which EC flux towers were maintained in the three post-clearcut chronosequences (BC-HDF, SK-HJP and QC-HBS) are described in Table 1 and in references cited therein. EC measurement methodology followed standard CCP protocol at each site, as described in these references and online at http://www.fluxnet-canada.ca/home.php?page= home&setLang=en. All towers had a fetch longer than 1 km, with footprints within the fetch except under very stable nighttime conditions when measurements were rejected. Before evaluating model performance against CO 2 fluxes measured at these sites, the model had first to reproduce site conditions by simulating site histories. This was accomplished by initializing ecosys at each site with the biological properties of the key tree and understory species, represented as plant functional types, and with the physical and chemical properties of the dominant soil type at each site (Table 1). Each plant functional type (coniferous or deciduous, tree or understory) had identical properties at each site in which it was present, except for those properties Table 1. Key attributes of vegetation, soils and climates at CCP post-clearcut chronosequences and reference forest sites   Tellus 62B (2010), 5 associated with timing of key phenological events such as leaf flush and reproductive growth (e.g. Grant et al., 2007a) which were adapted to the climate zone (temperate or boreal) at each site. All plant properties were derived from independent experiments and so remained unchanged from those used in earlier studies (e.g. Grant et al., 2007bGrant et al., , 2008Grant et al., , 2009a. Ecosys was also initialized with stocks of WD and fine litter estimated to remain after a stand-replacing fire in a model year corresponding to the late 18th or early 19th century, depending on site history. The model was then run from forest seeding to maturity (ca. 100 yr) under repeating sequences of continuous hourly weather data (radiation, air temperature, dewpoint or RH, wind speed and precipitation) recorded during the period of measurement at each site ( Table 1). The modelled forests were then clearcut or burned, depending on site history, in the model year corresponding to that in which the most mature forest stands (BC-DF49, SK-SOJP and QC-EOBS) were disturbed (Table 1). These disturbances used removal coefficients for foliage, other non-woody phytomass, living and standing dead wood derived from Kurz et al. (2009). All remaining phytomass was transferred to standing dead stocks, or to WD and fine non-woody litter stocks on the ground surface and in the soil profile (dead fine and woody roots). Standing dead stocks were transferred to surface WD stocks using a first-order hourly rate constant that gave a mean residence time in standing dead stocks of about 10 yr. Overstory and understory species were then reseeded and regrown under continuing sequences of weather data until they reached the age of the mature forest sites during the period of EC measurements.
To simulate the recently clearcut sites (BC-HDF00, BC-HDF88, SK-HJP02, SK-HJP94, SK-HJP75 and QC-HBS00), alternative runs were started at the beginning of the model year corresponding to that of the clearcut at each site (Table 1), using values for all state and driver variables stored at the end of the previous year from the model runs for the mature forest sites. The stands were then clearcut on harvesting dates reported from each site using removal coefficients for foliage, other nonwoody phytomass, living and standing dead wood derived from Kurz et al. (2009) or from local observations (e.g. Grant et al., 2007b) (Table 1). The stands were then reseeded with tree and understory species on planting dates reported from each site, and regrown under continuing sequences of hourly weather data until they reached the age of the clearcut forest stands during the period of EC measurements.
The model was thus run at all sites for ca. 200 yr before comparison with measured data, so that model results were independent of initial conditions. A background mortality rate of 1.2% (QC-EOBS) or 0.8% (all other sites) per year was applied to forest stands during the model runs, simulating natural self-thinning (Aakala et al., 2007). These rates were applied monthly to populations with trees of uniform biomass, so that both biomass and population were reduced to the same extent by each mortality event. All biomass removed by mortality was added to surface or subsurface detritus stocks in the model. During the last 150 yr of the model runs, atmospheric CO 2 concentration (C a ) rose exponentially from 280 to 385 μmol mol −1 , and precipitation NH 4 + and NO 3 − concentrations used to simulate wet N deposition rose exponentially from historical values based on Holland et al. (1999) to current values based on Meteorological Service of Canada (2004). Atmospheric concentration of NH 3 used to simulate dry N deposition was maintained at 0.0025 μmol mol −1 .

Model testing
Yearly regressions were conducted of modelled hourly CO 2 fluxes on EC hourly averaged CO 2 fluxes at each site in the three chronosequences for all years of measurement from 1 January 2001 to 31 December 2007. Model performance was evaluated from regression intercepts (a → 0), slopes (b → 1) and correlation coefficients (R 2 → 1), and from comparisons of root mean squares for differences between EC and modelled fluxes (RMSD) with root mean squares for error in EC fluxes (RMSE).
Values of RMSE at each site-year were calculated as the pooled root mean square of uncertainty in 1/2-hourly EC fluxes during the year using equations for CO 2 random flux measurement errors derived over comparable forest stands by Richardson et al. (2006). Values of NEP from ecosys were then compared with those derived from EC fluxes at hourly, daily and annual time scales, gap-filled according to established FCRN protocol described in Barr et al. (2004).
Comparisons were also made of modelled and measured C stocks in live aboveground tree biomass (bole, branches and foliage), live understory, and down WD (coarse, medium and small size classes >1 cm, but excluding stumps, standing dead, fine WD <1cm and surface LFH layers) in groundplots located in each site, which had been measured following guidelines established for Canada's National Forest Inventory (NFI Taskforce 2008). The simple average of multiple groundplots within each site was calculated for comparison and therefore may not reflect the spatial variation in C stocks within the tower footprint and true site mean as noted by Chen et al. (2009). Modelled near-surface soil temperatures were also compared with values measured at differently aged sites.

Model projections
Model sensitivity to changes in harvesting practices on ecosystem C stocks during subsequent forest regrowth was examined by running ecosys at BC-HDF88 for 80 yr after clearcutting in 1988 under repeating sequences of 1988-2007 weather data. Runs were conducted with removal coefficients for foliage, other non-woody phytomass, and living wood used in the model tests (Table 1), and then raised by 0.2 and 0.4 to simulate the effects on NEP of greater bole wood and aboveground detritus removal as might occur with whole tree and bioenergy harvests.

Tests of modelled versus measured CO 2 fluxes
Regression parameters of hourly CO 2 fluxes modelled versus measured during the entire year were poor (b < 0.8, R 2 < 0.5) for the recent clearcut sites during most years (BC-HDF00, SK-HJP02, SK-HJP94 and QC-HBS00 in Table 2). These parameters improved somewhat (1.0 > b > 0.8, R 2 > 0.6) in the 10-20year-old post-clearcut sites (BC-HDF88 and SK-HJP75), and further (1.2 > b > 1.0, R 2 > 0.75) in the older forest sites (BC-DF49, QC-EOBS). For all site-years, RMSD between modelled and measured hourly CO 2 fluxes was similar to or smaller than RMSE of EC fluxes calculated with the algorithm of Richardson et al. (2006) (Table 2). This indicated that much of the variance in EC fluxes not explained by the model could be attributed to random errors in flux measurements. Therefore better agreement between modelled and measured fluxes than that shown in Table 2 could not conclusively be attributed to improved model performance.
The poorer model performance in the more recently clearcut sites was attributed to low diurnal variation in EC CO 2 fluxes. Consequently variation in EC fluxes associated with surface boundary conditions and therefore amenable to modelling remained small relative to EC random error. Diurnal variation in EC fluxes was greater with respect to random error in the older forest sites (SK-HJP75, QC-EOBS and particularly BC-DF49), so that the model was able to explain a larger fraction of variation in the measured fluxes (Table 2).

Changes in hourly and daily CO 2 exchange with time after clearcutting
Both influxes and effluxes of hourly CO 2 measured and modelled during the period of most rapid CO 2 exchange in late June increased with time since clearcutting in the BC-HDF (Fig. 1), SK-HJP (Fig. 2) and QC-HBS (Fig. 3) chronosequences. In each of the recently clearcut sites (BC-HDF00 in Fig. 1, SK-HJP02 in Fig. 2 and QC-HBS in Fig. 3), CO 2 effluxes measured and modelled during the first 5 yr after clearcutting were driven by R h which increased gradually with colonization (eq. 10) and consequent accelerated decomposition (eq. 1) of the substantial amounts of WD and fine litter left from clearcutting. Decomposition at these sites was hastened by higher soil temperatures (Fig. 4) that arose in the model from solving the surface energy balance under increased exposure to incoming radiation caused by low tree and understory LAI. CO 2 effluxes were larger in the 10-20-year-old post-clearcut sites (BC-HDF88 in Fig. 1, SK-HJP94 in Fig. 2) and rose little further in the older forest sites (BC-DF49 in Fig. 1, SK-HJP75 in Fig. 2, QC-EOBS in Fig. 3), as declines in R h modelled from declining detritus stocks (eq. 10) and temperatures (Fig. 4) offset rises in R a modelled with tree growth (eq. 11).
CO 2 influxes rose gradually after clearcutting and reseeding from near zero values in the recently clearcut boreal sites (SK-HJP02 in Fig. 2 and QC-HBS in Fig. 3), through intermediate values in the 10-20-year-old post-clearcut sites (BC-HDF88 in Fig. 1, SK-HJP94 in Fig. 2), to maximum values in the older forest sites (BC-DF49 in Fig. 1, SK-HJP75 in Fig. 2 and QC-EOBS in Fig. 3). Influxes modelled during the first decade after clearcutting were mostly from pioneer deciduous plant functional types, much of which contributed to CO 2 effluxes through litterfall, while those modelled thereafter were mostly from dominant tree functional types (Table 1). In the model, these influxes rose because forest regrowth drove root elongation (eq. 16) and leaf expansion (eq. 20) that gradually raised nutrient uptake (eq. 17) and CO 2 fixation (eqs C1 to C11 from Grant et al., 2007b) through further elongation and expansion. Larger CO 2 influxes raised σ C and hence R a , thereby raising the autotrophic contribution to CO 2 effluxes. Influxes eventually stabilized in older forest stands as further root elongation and leaf expansion gave declining increments in nutrient uptake and CO 2 fixation but rising increments in senescence and litterfall (eqs 13, 15 and 18). Consequently declines in R h and rises in R a eventually ended with rising litterfall, mortality, and slowing growth, causing R e to stabilize in the older forest stands.
Modelled rises in CO 2 influxes were constrained in all three chronosequences by low leaf N concentrations (eqs A7e-k in Grant et al., 2010), calculated in the model as (M N +σ N )/(M C +σ C ) (eq. 19) ( Table 3). These concentrations were smaller than 0.029 g N g C −1 below which growth of Douglas fir was found to be N-limited by Hopmans and Chappell (1994). The early rise and subsequent decline in leaf N concentrations modelled with time since clearcut in the BC chronosequence was consistent with that measured in Douglasfir clearcuts by Bradley et al. (2002) and in Mountain hemlock by Titus et al. (2006). However the decline measured in the SK chronosequence was smaller than that modelled, possibly because of variation in soil quality among the chronosequence sites (table 1 in Mkhabela et al., 2009).
These low leaf N concentrations were caused in the model by competition for nutrient uptake with heterotrophic microbial populations colonizing WD and soil litter layers (eq. 10) with large C:N and C:P ratios. These ratios (Table 1) were larger than 28:1, above which growth of Douglas fir was found to become N-limited by Edmonds and Hsiang (1987). These ratios slowed decomposition and assimilation of N and P with respect to C (eqs. 2 and 6), forcing immobilization of mineral N and P (eq. 8b) and thereby slowing root and mycorrhizal uptake (eq. 17) and transfer to leaves (eq. A7b in Grant et al., 2010).

Site
Year  BC-HDF00 (Fig. 5). CO 2 influxes measured and modelled in recently clearcut sites generally failed to offset CO 2 effluxes (BC-HDF00 in Fig. 1, SK-HJP02 in Fig. 2 and QC-HBS00 in Fig. 3), causing these sites to be net sources of C during 2003-2007 (BC-HDF00 in Fig. 5, SK-HJP02 in Fig. 6 and QC-HBS00 in Fig. 7). CO 2 influxes measured and modelled in the 10-20-yearold post-clearcut sites more fully offset CO 2 effluxes at the two sites for which measurements were available (BC-HDF88 in Fig. 1, SK-HJP94 in Fig. 2), causing these sites to be nearly C neutral during this period (BC-HDF88 in Fig. 5, SK-HJP94 in Fig. 6). CO 2 influxes measured and modelled in the older forests were greater than CO 2 effluxes (BC-DF49 in Fig. 1, SK-HJP75 in Fig. 2 and QC-EOBS in Fig. 3), so that all older sites were net C sinks (BC-DF49 in Fig. 5, SK-HJP75 in Fig. 6 and QC-EOBS in Fig. 7). The seasonal pattern of modelled NEP also changed with time since clearcutting as the early pioneer vegetation dominated by deciduous functional types (BC-HDF00 in Fig. 5, SK-HJP02 and SK-HJP94 in Fig. 6 and QC-HBS00 in Fig. 7) were succeeded by sapling and then closed forests dominated by tree functional types (BC-HDF88 and BC-DF49 in Fig. 5, SK-HJP75 in Fig. 6 and QC-EOBS in Fig. 7). CO 2 fluxes modelled and measured during forest regrowth were also affected by weather. For example at the BC sites, influxes rose and effluxes declined with cooling during DOY 171-175 in 2004 while influxes declined and effluxes rose with warming during the same DOY in 2006 (Fig. 1). Consequently NEP at the BC sites was adversely affected by periods of warm weather, notably during the very warm summer of 2004 (Fig. 5). These effects of weather were greater in the older forest site (DF49) due to greater phytomass and greater hydraulic constraints in taller trees as described in Grant et al. (2007b).

Changes in annual respiration and productivity with time after clearcutting
Colonization (eq. 10) and decomposition (eq. 1) of WD and fine litter caused annual rates of modelled R h and R e to rise gradually for the first 5 yr after clearcutting (BC-HDF00, SK-HJP02 and QC-HBS00 in Table 4), to maintain larger values during the second decade after clearcutting (BC-HDF88 and SK-HJP94 in Table 4), and to rise little thereafter (BC-DF49, SK-HJP75, SK-SOJP and QC-EOBS in Table 4), consistent with EC-derived R h (Table 4) and CO 2 effluxes (Figs 1-3). In the model, the small rises in R e from the older stands were attributed to rises in R a , offset in some cases by declines in R h (e.g. BC-DF49 in Table 4). WD decomposition rates modelled at BC-HDF88 and SK-HJP94 remained low for several years after clearcutting while WD colonization was starting, but rose gradually to 0.05 and 0.025 yr −1 , respectively after about 10 yr as colonization progressed (Fig. 8), comparable to rates for wood derived in similar ecozones from litterbags (Trofymow et al., 2002) and for WD in other chronosequence studies (Naesset, 1999;Janisch and Harmon, 2002;Shorohova et al., 2008;Melin et al., 2009). These rates allowed WD modelled at BC-HDF88 to approach full colonization after 20 yr (colonized WD approaching total remaining WD in Fig. 8), consistent with the time course of WD decomposition stages estimated in Ranius et al. (2003). Most WD at BC-HDF88 was observed to have reached decomposition class 3 in 2002 (Trofymow, FCRN DIS at http://fluxnet.ccrp.ec.gc.ca), consistent with the partial colonization modelled at this site in 2002 (Fig. 8).
Modelled and EC-derived GPP were much smaller than R e in recently clearcut sites (BC-HDF00, SK-HJP02 and QC-HBS00 in Table 4), but rose more rapidly with time, approaching R e in the 10-20-year-old post-clearcut sites (BC-HDF88 and SK-HJP94 in Table 4), and surpassing R e in the older forest stands 30-105 yr after disturbance (BC-DF49, SK-HJP75, SK-SOJP and QC-EOBS in Table 4). The greater rise in GPP versus R e with time since clearcutting arose from the greater rises in CO 2 influxes versus effluxes (Figs 1-3), causing rises in daily NEP (Figs. 5-7) that caused those in annual NEP (Table 4). These rises were accelerated by fertilization of the BC chronosequence in 2007, although EC-derived NEP rose less than did modelled NEP because fertilization raised EC-derived R e at BC-HDF00 and BC-HDF88 and lowered EC-derived GPP at BC-DF49 (Table 4) even though 2007 was a comparatively cool year at these sites (MAT = 7.7 • C).  stocks that were compared with similar aboveground C stocks measured in groundplots at each site (Table 5). During model runs, surface WD stocks (excluding standing dead) were replenished by detritus from harvest removals (Table 1), and by ongoing litterfall from living and standing dead phytomass (eq. 13), and depleted by R h (eq. 3; Table 4). The particularly large WD stock recorded at BC-HDF88 in 2004 was modelled by using a comparatively small harvest removal co-efficient at this site in 1988 (Table 1). Modelled WD stocks tended to exceed the site average groundplot values for all sites in SK and QC (Table 5), Much of the detritus modelled after clearcutting was carried forward from the previous forests (e.g. Fig. 8) and thus in an advanced stage of decomposition, and so may not have been adequately measured in the groundplots. Differences in the definitions of modelled and measured pools could also account for some of the discrepancies.

Changes in ecosystem C stocks with time after clearcutting
Plant growth driven by GPP (Table 4) generated live tree (bole, branch and foliage) and understory (foliage and branch) phytomasses during model runs that increased slowly during the first and second decades after clearcutting, but more rapidly thereafter, consistent with the slow but continuous rise in CO 2 influxes modelled and measured with time since clearcutting (Figs 1-3). These phytomasses were, with some exceptions, usually within one standard deviation of the average of measured groundplot values at each site (Table 5), even with common model parameter values used at all sites for all biological properties of each plant functional type. Chronosequences represent a substitution of space for time and although all sites in a chronsequence were modelled as having the same site conditions and pre-clearcut histories, in reality these differed among sites which could further account for the discrepancies between modelled and measured values. Tree wood growth modelled over 100 yr after clearcutting tracked growth curves derived from inventory data for the ecozone in which each chronosequence was located ( fig. 12 in Grant et al., 2009a).

Changes in net ecosystem productivity with changes in harvest removal
The effects of fine litter and WD stocks left after clearcutting on the subsequent time courses of R h , GPP and hence NEP suggest that the harvest practices that determine these stocks might affect subsequent forest productivity. In the model, these stocks depended on forest age and productivity before clearcutting, and also on removal coefficients for foliage, other nonwoody phytomass, living wood and standing dead wood used to represent the effects of logging practices on forest C distribution (Table 1). The modelled responses of R h , GPP and NEP to changes in removal coefficients were examined by projecting forest growth for 80 years after clearcutting at BC-HDF88 with removal coefficients for foliage, other non-woody phytomass, and live wood raised by 0.2 and 0.4 from those used in the model tests.
Raising these coefficients reduced the amounts of WD and fine litter left in the model after clearcutting, thereby reducing initial soil + detrital C stocks (Fig. 9a). However these reductions caused slower R h thereafter (eq. 3), so that soil + detritus C stocks gradually converged towards common values after 80 yr. During the first several years after clearcutting, slower R h from greater harvest removal caused less immobilization (eq. (8)) of N mineralized during detrital decomposition (eq. 2a), thereby raising N availability for root uptake. However slower R h from greater harvest removal also caused less asymbiotic N fixation (eqs A26, A27 in Grant et al., 2007b), causing less mineral-ization of diazotrophic decomposition products, which reduced the gains in N availability from slower immobilization. Raising harvest removal coefficients also reduced the thickness of the surface WD layer, thereby raising temperatures modelled in the forest floor underneath (Fig. 10), as found experimentally by Tan et al. (2009). This warming hastened N mineralization from decomposition (eq. 2a), further raising N availability with greater harvest removal.
The time course of N uptake and assimilation caused foliar N contents in the model to rise during decomposition and mineralization of low C:N fine litter during the first five years after clearcutting (assart effect) ( Fig. 9b; Table 3). Foliar N contents then declined during slower decomposition and mineralization of high C:N WD over the following 25 yr, and then rose gradually thereafter during decomposition and mineralization of lower C:N products of litter and WD decomposition. Greater N availability from greater harvest removal hastened early plant N uptake in the model (eq. 17), causing foliar N concentrations to rise slightly during the first 30 years after clearcutting (Fig. 9b). However greater harvest removal caused an eventual reduction in mineralization from decomposition products which caused foliar N concentrations to decline slightly thereafter. These changes in foliar N concentrations with harvest removals slightly hastened tree phytomass growth for several decades after clearcutting, but slowed it thereafter, so that net effects of harvest removal on tree phytomass were small (Fig. 9c). Some of the early gain in phytomass modelled with greater harvest removal occurred in the deciduous understory species rather than in the trees. However the modelled understory eventually declined with rising tree LAI, so that most of its C and N stocks were remineralized by 20 yr after clearcutting in all the harvesting treatments. The net effect of greater harvest removal on forest growth in the model thus arose from complex interactions among N transformation processes on N availability that changed over time.
These model results were consistent with findings from other studies that the intensity of wood and litter removal with clearcutting have variable effects on subsequent tree growth. Tan et al. (2009) found that whole tree harvest with slash removal stimulated early lodgepole pine and Douglas-fir growth compared to stem-only harvesting at some BC sites as modelled in our study, but not at others. Olsson and Staaf (1995) found that removal of dense slash layers left from clearcutting boreal pine hastened the early growth of some plant species but not others, although this effect declined as slash decomposed over time, as also modelled in our study. Holcomb (1996) observed that removal of logging slash and forest floors after harvest may reduce productivity at sites with little SOC accumulation, but raise it at sites with greater accumulations, such as modelled at BC-HDF88 in our study. Runs of the FORCYTE-11 process model, tested and calibrated for a Douglas-fir site in the Shawnigan Experimental forest on southern Vancouver Island, predicted that with whole tree harvest, wood production would decline 10-30% after multiple rotations, especially if rotations were short, as accelerated N export of N rich foliage and juvenile wood was greater than the annual N deposition for this site (Trofymow and Sachs, 1991). These variable results in the literature and in our modelling study are consistent with observations that the intensity of wood and litter removal after clearcutting have little effect on C and N pools (Olsson et al., 1996) or on N mineralization rates (Titus et al., 2006) in the soil organic and mineral layers from which most plant N uptake occurs.
By slowing R h (Fig. 9a) and hastening phytomass growth (Fig. 9c), greater harvest removals caused earlier rises in Table 4. Annual net ecosystem productivity (NEP), gross primary productivity (GPP), ecosystem respiration (R e ) and heterotrophic respiration (R h ) derived from gap-filled eddy covariance measurements (EC)  Year  Fig. 8. Colonized and total woody debris (WD) (= S C and S C + S C in eq. 10) remaining in the model to 2007 after clearcutting at BC-HDF88 in 1988 andSK-HJP94 in 1994. modelled NEP, allowing C neutrality to be reached several years sooner (Fig. 9d). This model response was consistent with that of an empirical model based on inventories in post-disturbance chronosequences of Douglas fir-dominated forests by Janisch and Harmon (2002). In their model, greater wood removal reduced the early decline in NEP and hastened the subsequent rise in NEP to C neutrality after clearcutting. However NEP in ecosys converged to similar values during later forest regeneration more than 30 yr after the different harvest removals (Fig. 9d), in contrast to NEP predicted during regeneration by Janisch and Harmon (2002) which rose less and remained lower with greater removal. This convergence in ecosys occurred because removal effects on C stocks in soil + detritus (Fig. 9a) and trees (Fig. 9c) both diminished over time, and so only affected NEP during the first three decades after clearcutting. Model results collectively indicated that greater harvest removal at BC-HDF88 would have a large effect on NEP within 30 yr after clearcutting (Fig. 9d), but only a limited effect thereafter. Raising removal coefficients by 0.2 and 0.4 from those used in model testing caused harvest removals to increase by 3645 and 7290 g C m −2 , respectively. These increased removals caused reductions in ecosystem C stocks (living + standing dead trees, WD, fine litter, humus, and losses through DOC + DIC export) of 541 and 1188 g C m −2 , respectively after an 80-yr harvest cycle, or about 15% of the increased removals. These losses represented the cost of the removals to ecosystem C stocks at the following harvest.

Controls on the time courses of respiration and productivity after clearcutting
The time courses of NEP modelled and measured after clearcutting at the three chronosequences in this study arose from different time courses of respiration and primary productivity. Modelled and EC-derived R h rose with higher detritus stocks (Table 5) and temperatures (e.g. Fig. 4) over 2-7 yr of measurement in the recently clearcut sites at BC-HDF00, SK-HJP02 and QC-HBS00, remained high in the 10-20-year-old post-clearcut sites at BC-HDF88 and SK-HJP94, then declined and eventually stabilized after more than 20 yr with declining detritus stocks and temperatures and rising litterfall and mortality in the older forest sites at BC-DF49 and SK-HJP75 (Table 4). Some interannual variability was superimposed on this trend, such as higher values modelled and measured during 2004, a year with a particularly warm summer at the BC chronosequence, and from 2004 to 2006, during which warming occurred at the SK and QC chronosequences (Grant et al., 2009a). This time course of R h with time since clearcutting in the model appeared BC-HDF00 (2002) Martin et al. (2002) found that decomposition rates measured in a post-clearcut silver firwestern hemlock chronosequence in BC rose during the first 6 yr, and declined thereafter. Kolari et al. (2004) found that soil respiration R s measured in a post-clearcut pine chronosequence in southern Finland was highest in a 12-year-old stand and lower in both younger and older stands. However in some cases drier surface conditions were found to slow decomposition at exposed clearcut sites in post-harvest chronosequences of Douglas-fir on southern Vancouver Island (Trofymow, 1998;Addison et al., 2003). Declining contributions of R h to R e more than 5 yr after clearcutting may therefore be a general phenomenon, causing R e to rise little with forest age more than 20 yr after clearcutting (Table 4). In the model, the time course of R h was simulated with a function for microbial colonization of WD and fine litter (eq. 10) that simulated gradual rises in decomposition and hence R h from new detritus for several years after clearcutting followed by gradual declines (Fig. 8), as frequently observed in field studies (e.g. Mäkinen et al., 2006;Montes and Cañellas, 2006). The time course of R h in the model caused NEP to remain negative but stable for the first decade after clearcutting, and to rise only gradually during the second decade. This time course of NEP generally matched that of the EC-derived NEP (Figs. 1-3, 5-7; Table 4). However this time course could not be modelled from simple first-order functions of the amounts of WD and fine litter left after clearcutting as done in many empirical models, which can give only declining R h and hence rising NEP with time following clearcutting.
In contrast to the time course of R e as driven by R h , that of GPP indicated a gradual but continuous rise from reseeding to maturity (Figs 1-3) which drove the gradual rise in NEP beginning several years after clearcutting (Table 4). The rate at which GPP rose in the model was constrained by low foliar N concentrations (Table 3), the time course of which was similar to that measured in other field studies. Bradley et al. (2002) found that foliar N concentrations declined from 24.4 to 17.6 and 16.3 mg N g C −1 (assuming needles are 50% C) in 4-, 7-and 11year-old stands, respectively in a post-clearcut chronosequence of silver fir-western hemlock in coastal BC. Bradley et al. (2002) attributed this decline to measurements of mineral N availability that rose to maximum values 4-5 yr after clearcutting, and declined thereafter. This time course of mineral N availability was found to follow that of slash decomposition after clearcutting as measured by Martin et al. (2002) and modelled in this study (Table 4).
The time course of foliar N concentrations was hypothesized in the model to arise from the time courses of decomposition (eqs 1, 2; Fig. 8), microbial growth (eqs 7, 9), colonization (eq. 10), nutrient mineralization versus immobilization (eq. 8) and asymbiotic N 2 fixation (eqs A26, A27 in Grant et al., 2007b) in WD and fine litter complexes with, respectively higher versus lower C:N ratios, and lower versus higher specific decomposition rates. These processes controlled the time course of soil  mineral N concentrations that drove root and mycorrhizal N uptake (eq. 17) and plant N assimilation (eq. 19; Table 3) that in turn drove growth in root and mycorrhizal lengths (eq. 16) and in leaf areas (eq. 20), and hence rises in N uptake (eq. 17) and GPP (eq. A7 in Grant et al., 2010).
Interactions among these processes generated time courses for long-term projections of C accumulation in pioneer shrubs and dominant trees that reproduced the Chapman-Richards functions used to estimate slow early growth in empirical models (e.g. Janisch and Harmon, 2002;Taylor et al., 2005;Kurz et al., 2009) under diverse site conditions (Table 5; Grant et al., 2009a). When the time courses of GPP that drove C accumulation were combined with those of R e as driven by R h , the model was able to generate different time courses of NEP among the three post-clearcut chronosequences that were consistent with EC-derived values, without resorting to site-and speciesspecific parametrizations. Thus the more productive temperate BC sites with greater WD stocks left by clearcutting (Table 5) had greater R h and lower NEP for several years after clearcutting, but attained greater GPP and hence NEP during later years than did the less productive boreal sites in SK and QC (Table 4).
Continuous EC measurements of CO 2 fluxes along diverse post-clearcut chronosequences have improved constraints to testing a process-based model of decomposition and regrowth. This testing has enabled a robust simulation of post-harvest changes in forest C stocks, allowing the model to be used in lifecycle studies of forest C following harvest (e.g. Fig. 9) without the need for site-specific parameterization.

Acknowledgments
Computing facilities for ecosys were provided by the University of Alberta and by Cybera, a corporation managing cyberinfrastructure-related technologies in collaboration with academic and industry partners. Funding for this study was provided by a CFCAS grant to the Canadian Carbon Program (CCP), National Sciences and Engineering Research Council of Canada (NSERC) network and strategic grants to the University of British Columbia (UBC) and by the Natural Resource Canada PERD program. We thank Bob Ferris, Gurp Thandi, Colin Ferster, Mark Gillis and Frank Eichel of the Canadian Forest Service (CFS) and staff of B.A. Blackwell and Associates for their help collecting and processing National Forest Inventory-style ground plot data for the BC coastal sites. g C g C −1 h −1 (1) 1.0 (protein), 1.0 (carbo.), 0.15 (cellulose), 0.025 (lignin) Grant et al. (1993a,b) E n Energy requirement for growth of M C kJ g C −1 (5) 25 F Fraction of microbial growth allocated to kinetic components of M C -(7, 9) 0.55 (labile), 0.45 (resistant) Grant et al. (1993a,b) K iD Inhibition constant for M C on S C g C m −3 (1) 25 Grant et al. (1993a,b Grant and Hesketh, 1992) β Specific colonization of detritus by growth of M C g C g C −1 h −1 (10) 2.0 θ P Root or mycorrhizal porosity m3 m-3 (16) 0.25 Grant (1998)