Climate and atmospheric drivers of historical terrestrial carbon uptake in the province of British Columbia , Canada

The impacts of climate change and increasing atmospheric CO2 concentration on the terrestrial uptake of carbon dioxide since 1860 in the Canadian province of British Columbia are estimated using the process-based Canadian Terrestrial Ecosystem Model (CTEM). Model simulations show that these two factors yield an enhanced carbon uptake of around 44 gC m−2 yr−1 (or equivalently 63 gC m−2 yr−1 over the province’s forested area), during the 1980s and 1990s, and continuing into the 2000s. About three-quarters of the simulated sink enhancement in our study compared to pre-industrial conditions is attributed to changing climate, and the rest is attributed to increase in CO 2 concentration. The model response to changing climate and increasing CO 2 is corroborated by comparing simulated stem wood growth rates with ground-based measurements from inventory plots in coastal British Columbia. The simulated sink is not an estimate of the net carbon balance because the effects of harvesting, insect disturbances and land-use change are not considered.


Introduction
Atmospheric carbon dioxide (CO 2 ) concentration is increasing due to emissions from anthropogenic use of fossil fuels and changes in land use.About half of the emitted anthropogenic carbon is taken up by land and ocean, and this uptake has slowed the rate of increase of atmospheric CO 2 in response to anthropogenic emissions (Canadell et al., 2007;Le Quéré et al., 2013).Over land, several lines of evidence indicate that carbon uptake is occurring in northern mid-to high-latitude regions (Ciais et al., 2010).This evidence includes (1) flux towers measuring the land-atmosphere exchange of CO 2 typically over an area of a few square kilometres (e.g.Krishnan et al., 2008;Yuan et al., 2009), (2) inversion-based studies using atmospheric transport models in conjunction with observations of atmospheric CO 2 to infer the location of sinks and sources of carbon on the continental scale (e.g.Deng et al., 2007;Deng and Chen, 2011;Gourdji et al., 2012), (3) inventory-based studies (Pan et al., 2011), and (4) modelling approaches where terrestrial ecosystem models are driven by observed atmospheric CO 2 concentration and climate (e.g.Huntzinger et al., 2012).A combination of forest inventory-based assessments and terrestrial ecosystem modelling has been used for assessing carbon uptake over land in Europe (Luyssaert et al., 2010) and Canada (Stinson et al., 2011).
The various approaches have their characteristic strengths and weaknesses.For example, the spatial coverage of flux towers is too small to infer regional-or continental-scale land-atmosphere CO 2 exchange, although when used in conjunction with meteorological data their measurement can be extrapolated to the global scale (e.g.Beer et al., 2010).In inversion-based studies the inferred carbon sink is a function of atmospheric transport of CO 2 and prior estimates of CO 2 fluxes and boundary conditions (see e.g.Gourdji et al., Published by Copernicus Publications on behalf of the European Geosciences Union.2012).However, the inferred fluxes do provide useful information at continental scales.The inventory-based estimates (e.g.Pan et al., 2011) are able to provide large-scale estimates of carbon sink and source distribution and are able to account for the age-class structure of forested landscapes and the strong, age-dependent impacts on net ecosystem production, but there are large differences in available data, with data limitations especially for tropical regions.Finally, modelling approaches depend on assumptions made in the models, and, as a result, estimates vary from model to model even when driven by identical forcing.Huntzinger et al. (2012), for example, find that the simulated net atmosphere-land CO 2 exchange over North America for the 2000-2005 period ranges from −0.7 Pg C yr −1 (a source) to 2.2 Pg C yr −1 (a sink) across nineteen terrestrial ecosystem models.Regionalscale modelling studies fill in the spatial gap between the point scale flux tower measurements and large continentalscale inversion-based studies.Models typically provide results at higher spatial resolution than inversion-based studies.In the end, these various approaches are complementary, and together provide knowledge that helps reduce the overall uncertainties in regional-to global-scale estimates of terrestrial carbon exchange.
Here we use the process-based Canadian Terrestrial Ecosystem Model (CTEM) coupled to the Canadian Land Surface Scheme (CLASS) to simulate the terrestrial carbon budget over the province of British Columbia (BC) in Canada from 1900 to 2010.Our primary objective is to estimate the magnitude of the terrestrial carbon sink response to changes in climate and increasing atmospheric CO 2 concentration, which stimulates plant growth through increased photosynthesis rates (i.e. the CO 2 fertilization effect).To our knowledge this is the first such assessment for BC using a processbased terrestrial ecosystem model, although such models have been employed in the past at the global (Ahlström et al., 2013), continental (e.g.Huntzinger et al., 2012;Zhang et al., 2012) and regional scales (e.g.McGuire et al., 2010).Located between 48 • and 60 • N, the terrestrial ecosystems of the Canadian province of BC are representative of mid-to high-latitude coniferous forests and represent a wide range of conditions from coastal to interior climates, from very wet to very dry conditions, and from low-to high-elevation ecosystems.The insight gained into response of terrestrial ecosystems in the province to climate and CO 2 forcing is likely to apply to other mid-to high-latitude forests.
Estimates of the carbon budget in BC forests from 1920 to 1989 have been previously obtained using version 2 of the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS2) (Kurz et al., 1996), while estimates for managed forests, including those in BC, from 1990 to 2008 have been obtained using an updated version (CBM-CFS3; Kurz et al., 2009) of this inventory-based model (Stinson et al., 2011).Tree growth in the CBM family of models, however, only implicitly accounts for changes due to the CO 2 fertilization effect and the impacts of climate change which are inferred empirically from past tree growth (Hember et al., 2012) as reflected in empirical yield curves.The processbased model used here, on the other hand, simulates ecosystem response of all plant functional types (PFTs), including trees, and explicitly accounts for changes in ecosystem processes due to CO 2 fertilization and climate change.
In Sect. 2 we describe the modelling framework used in this study, including the model subcomponents and the forcing data sets.Our main results are presented in Sect.3. Finally, our results are discussed more broadly in Sect. 4 and a summary of our main conclusions is provided in Sect. 5.

Process-based terrestrial ecosystem model
The process-based terrestrial ecosystem model used in our study is CTEM (Arora and Boer, 2005a, b), the interactive vegetation component of the Canadian Centre for Climate Modelling and Analysis (CCCma) Earth System Model (CanESM2; Arora et al., 2011), which is coupled to version 3.5 of CLASS, the physical land surface component of CanESM2 (Verseghy, 1991;Verseghy et al., 1993).The coupled system is driven by half-hourly meteorological data to simulate the physical (including soil moisture, soil temperature and snow) and biogeochemical (including vegetation biomass and soil carbon) states of the land surface.The primary terrestrial ecosystem processes that are modelled in CTEM for this study are photosynthesis; autotrophic and heterotrophic respiration; allocation; phenology; turnover of leaves, stems and roots; and fire.CTEM simulates vegetation growth and calculates time-varying carbon storage in three living vegetation pools (leaves, stems and roots) and two dead carbon pools (litter and soil organic matter).When coupled with CLASS, CTEM estimates canopy conductance which is then used for energy and water balance calculations in CLASS.A single-leaf photosynthesis approach is used with coupling between photosynthesis and canopy conductance based on vapour pressure deficit (Leuning, 1995).Photosynthesis and leaf maintenance respiration calculations are performed half-hourly, while other biogeochemical processes are simulated at a daily time step.While CTEM simulates fire, the effects of forest harvesting and insect disturbances are not modelled.
In CTEM's standard configuration, used in CanESM2 for global application, nine PFTs are considered.However, in this study, an additional needleleaf evergreen PFT is considered for simulating terrestrial ecosystem processes more realistically for the province of BC, and this is briefly discussed in Sect.2.4.The land surface scheme in CLASS employs the mosaic approach to account for sub-grid variation in vegetation types (Li and Arora, 2012), with each grid cell represented by a number of mosaic tiles depending on the PFTs present in a given grid cell.Biophysical energy, water balance and biogeochemical carbon balance calculations are performed over each mosaic tile, and the results are weighted by the areal fraction of each PFT to yield a grid cell mean value.There is no representation of age class or time since disturbance for the PFTs and all terrestrial ecosystem processes are modelled for an average-aged tree in the landscape.
The net atmosphere-to-land CO 2 flux F L (Tg C yr −1 ) in the CTEM is modelled as where H L = H V + H S is the total terrestrial carbon stock (Tg C) which is made up of living vegetation biomass in leaf, stem and root carbon pools (H V ) and dead organic matter carbon in soil and litter pools (H S ).Positive values of F L indicate that land is gaining carbon from the atmosphere.N is the terrestrial net primary productivity obtained as the difference between gross primary productivity (G) and autotrophic respiration (R A ), R H is the heterotrophic respiration and F F is the emissions associated with fire.The time-integrated version of Eq. ( 1) relates the change in total land carbon H L to the cumulative atmosphere-to-land CO 2 flux ( FL ) as Simulated F L is thus the net biome productivity (NBP) but without the carbon sources from harvest (F H ) and insect disturbance (F I ).Harvesting typically includes removal of timber, while insect disturbances increase the transfer of biomass to dead organic matter.Any residual dead organic matter after harvest or insect disturbance decomposes over time so these disturbances also affect R H .When estimates of F H and F I are available, net biome productivity F L is modified as (3)

Soil, vegetation and climate data
Soil depth, soil sand and clay content, and areal fraction of each PFT are based on data obtained from the Pacific Climate Impacts Consortium (PCIC), and are available at 1/16th-degree resolution (Schnorbus et al., 2011).Soil types were derived from the Global Soil Data Products (GSDT, 2000), which are built from the global pedon database produced by the International Soil Reference and Information Centre (ISRIC) (Batjes, 1995) and the FAO-UNESCO Digital Soil Map of the World (DSMW) (FAO, 1995).Three soil layers are considered, with the depth of the first and second layers uniform across the study area at 0.1 m and 0.25 m, respectively; the depth of the third layer varies spatially and depends on elevation and slope (Schnorbus et al., 2011).Fractional vegetation cover for each PFT in each grid cell was based on a 25 m resolution land cover data set from the Canadian Forest Service's Earth Observation for Sustainable Development of Forests (EOSD) (Wulder et al., 2003).Missing data that included no data or corruption by cloud and shadow (comprising 0.12 %, 0.50 % and 1.37 % of the land area, respectively) were filled using the University of Maryland's Advanced Very High Resolution Radiometer (AVHRR)-based global land cover classification data set (http://glcf.umd.edu/data/landcover/data.shtml).The EOSD data have been validated over Vancouver Island in BC with an accuracy of 86 % for coniferous forests (Wulder et al., 2008).The areas of forest classes mapped in the EOSD data set for year 2000 have also been corroborated by the areas in the 2001 Canadian Forest Inventory (Power and Gillis, 2006).Soil and vegetation data were interpolated to the spatial resolution of about 0.35 • latitude × 0.61 • longitude used in this study, which is equivalent to about 40 km × 40 km.The primary PFTs that exist in BC are needleleaf evergreen trees, broadleaf cold deciduous trees and C 3 grasses (see Table 1).A small fraction of the province is also covered with crops, but these were treated as C 3 grasses.Figure 1 shows the spatial distribution of the fractional vegetation cover of these three PFTs and the total fractional vegetation cover.At the ∼ 40 km resolution used here, some grid cells at the province's boundary lie partially outside its borders (see Fig. 1).As a result the area of the province used in the model is 1 005 388 km 2 , or about 6 % greater than the actual area of 944 700 km 2 .
Climate data to drive our model were obtained from the CRUNCEP data set , which is based on the National Centers for Environmental Prediction (NCEP) reanalysis (Kanamitsu et al., 2002) with monthly means adjusted to match the Climate Research Unit (CRU) observations.The CRUNCEP data (surface temperature, pressure, precipitation, wind, specific humidity, shortwave and longwave radiation fluxes) are available at a resolution of 0.5 degrees and at a six-hourly time interval (https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm4.CRUNCEP.v4.TPHWL6Hrly.html).Data were extracted for BC and spatially interpolated to ∼ 40 km × 40 km grid used in this study.The six-hourly data were disaggregated to half-hourly time resolution following the approach of Arora and Boer (2005a).Temperature, longwave radiation, wind speed, specific humidity, and pressure were linearly interpolated in time; shortwave radiation was assumed to change with solar zenith angle with a maximum value at the local solar noon; precipitation was disaggregated using the six-hourly precipitation amount to estimate the number of wet half hours, and total 6-hourly precipitation  amount was assumed to be randomly distributed over the wet half hours.
Figure 2  In particular, the weak positive trend in precipitation over BC is inconsistent with the analysis of station-based data by Mekis and Vincent (2011), who find that more than 80 % of the stations in the province of BC show trends of a greater than 10 % increase in precipitation over the period 1950-2009.Even larger trends are seen especially in the southern part of the province.The implications of this limitation in the CRUNCEP data are discussed later.

Observation-based data sets
Observation-based data sets of leaf area index (LAI), gross primary productivity (GPP), vegetation biomass and stem wood growth rates are used to assess modelled results.The observation-based estimates of LAI are based on remotely sensed data from Deng et al. (2006) (Ruesch and Gibbs, 2008).The second observation-based estimate of vegetation biomass is for commercial forests in BC (Penner et al., 1997) and reports an aboveground vegetation biomass density of 15.8 kg m −2 for managed forests, which we convert to carbon density (7.9 kg C m −2 ) by multiplying by 0.5 (Lamlom and Savidge, 2003;Sarmiento et al., 2005;Fonseca et al., 2011).The third observation-based estimate is derived from Canada's National Forest Carbon Monitoring, Accounting and Reporting System (NFCMARS), which gives an aboveground biomass of 5181 Tg C, or an average density of 7.8 kg m −2 (updated from Stinson et al., 2011), for all of BC's forests, averaged over the period 1990 to 2011.
Finally, we also use the observation-based stem wood growth rate from Hember et al. (2012), who analyse data from 1267 permanent inventory plots from throughout coastal BC between 1950 and 2002 that were installed by the British Columbia Ministry of Forests, Lands and Natural Resource Operations (BC MFLNRO).These measurements were made periodically over five-or ten-year intervals which allow for analysis of multi-annual and decadal variability.

Simulations
We used 1901 to 1940 meteorological data repeatedly and the 1860 CO 2 concentration value of 286.2 ppm to spin up and bring the model's carbon pools into equilibrium.We refer to this as our pre-industrial simulation.Transient historical simulations were then performed from 1861 to 2010 driven by atmospheric CO 2 increasing from 286.3 ppm in 1861 to 389.1 ppm in 2010.For the period 1861-1900 we used 1901 to 1940 climate to drive the CLASS and CTEM models, and for the period 1901-2010 we used the actual climate data.The fractional vegetation cover is kept constant over the duration of the simulation (1861-2010).The effect of land-use change is thus neglected but is not expected to affect our results significantly since farm land (comprising of both pasture and crop land area) covers only about 2.8 % of the province (Statistics Canada, 2011; http://www.statcan.gc.ca/pub/95-640-x/2012002/prov/59-eng.htm).Three transient historical simulations were performed that are driven by (1) increasing atmospheric CO 2 but with 1901-1940 climatology (the CO 2 simulation), (2) changing climate but with fixed pre-industrial CO 2 of 286.2 ppm (the CLIM simulation) and (3) increasing atmospheric CO 2 and changing climate (the CLIM+CO 2 simulation).These three simulations allow estimation of the separate contributions of CO 2 fertilization and climate change to the terrestrial uptake of carbon in BC.With a constant climate and atmospheric CO 2 concentration, terrestrial C stocks would remain unchanged from their equilibrium state.Therefore, any change away from the equilibrium carbon pools is attributable to changing atmospheric CO 2 concentration, changing climate or the combination of the two.In essence, the G, R A , R H and F F terms in Eq. ( 1) respond to changes in climate and atmospheric CO 2 concentration resulting in a positive or negative F L (i.e. a land carbon sink or a source) with associated changes in H V and H S pools.
Most of BC is covered with coniferous, needleleaf, evergreen trees.However, tree species in the interior region, which experience colder winters and drier summers than those in much of the coastal region, are known to be colder and more drought-tolerant compared to those in the coastal region.In its standard configuration used in CanESM2 for global application, CTEM considers nine PFTs (Table 1), which include one needleleaf evergreen PFT.However, for application over BC at a spatial resolution of 40 km we found that, while CTEM's single needleleaf evergreen PFT yields reasonable LAI and gross primary productivity (GPP) in the coastal region compared to observation-based estimates, it yields unrealistically low LAI and vegetation biomass in the interior of the province (see Appendix A).For the purposes of this modelling study, and as explained in the Appendix, we therefore make the first-order distinction based on climate and split CTEM's default needleleaf evergreen PFT into coastal and interior types.The coastal type retains the default parameter values of CTEM's needleleaf evergreen PFT.For the interior needleleaf evergreen PFT, we assume lower rates of leaf loss from cold and drought in the phenology parameterization of CTEM (Arora and Boer, 2005a) and assume a longer leaf life span following Reich et al. (1995) (see Appendix A).

Results
In Table 2 we compare modelled results with observationbased estimates.Modelled results are reported for a reference period which we choose to be 1990-2005 as well as for the time period that corresponds to each observation-based estimate.

LAI, GPP and vegetation biomass
Simulation results obtained by including the modified PFT set (which considers the interior and coastal needleleaf evergreen PFTs separately) are shown in Fig. 3 and Table 2, which compare simulated LAI and GPP with observationbased estimates from Deng et al. (2006) and Beer et al. (2010), respectively.The introduction of a second needleleaf evergreen PFT increases both LAI and GPP in the interior of the province and significantly improves comparison between simulated and observation-based LAI and GPP (compare Fig. A1 in the Appendix and Fig. 3).Limitations still remain, however, and in particular simulated LAI is lower in southern BC.This is because of either biases in the CRUNCEP data or the inability of CTEM to simulate vegetation realistically with just two needleleaf evergreen PFTs.A discussion of the challenge associated with parameterizing CTEM at the species level is presented in the Appendix.Nevertheless, the simulation with the modified PFT set yields a province-wide averaged GPP of 618 g C m −2 yr −1 (1998-2005 average) and province-wide areally averaged LAI of 2.4 m 2 m −2 (2000-2005 average) consistent with the observation-based estimates of 597 g C m −2 yr −1 for GPP (Beer et al., 2010) and of 2.5 m 2 m −2 for LAI (Deng et al., 2006), respectively (see Table 2).
Figure 4 shows the distribution of simulated living (vegetation) and dead (litter and soil carbon) pools.The spatial distribution of the vegetation biomass is similar to that for LAI, as expected.The spatial distribution of litter and soil carbon is somewhat different from that of vegetation biomass, with higher values towards the northern part of the province despite low vegetation density in the north.Colder climates yield lower litter and soil carbon decomposition rates and hence larger pool sizes.Table 2 also compares the simulated vegetation biomass with the three available observation-based estimates.The province-wide total vegetation biomass is 3901 Tg C averaged over the reference period of 1990-2005.The simulated forest-only biomass is 3714 Tg C which yields a forest-only biomass density of 5.6 kg C m −2 (when divided by BC forest area of 663 981 km 2 ) over the same time period.In Table 2, our estimates of total forest vegetation biomass and density are lower than the inventory-based estimates for BC from Penner et al. (1997) and Stinson et al. (2011), but our province-wide total vegetation biomass is almost twice that of the values obtained from the global data set of Ruesch and Gibbs (2008).Our simulated vegetation biomass thus lies between the estimates from Ruesch and Gibbs (2008) and the two inventory-derived estimates of Penner et al. (1997) and Stinson et al. (2011).

Net atmosphere-land CO 2 flux
Figure 5a shows the BC-averaged net atmosphere-land CO 2 flux F L over the 1900-2010 period for our three simulations (CO 2 , CLIM and CLIM+CO 2 ).Positive values indicate that land is gaining carbon from the atmosphere.Our results estimate that climate change and increasing atmospheric CO 2 have contributed a net sink of carbon since about 1930 and this sink grew to about 44 g C m −2 yr −1 during the 1980s, 1990s and early 2000s (the brown line in Fig. 5a from the CLIM+CO 2 simulation).In Fig. 5a, the temporal pattern of increasing F L that gradually levels off is the typical response of terrestrial ecosystems to a forcing that leads to carbon uptake (e.g.see Fig. 1c in Arora et al., 2013).As the vegetation and soil carbon mass slowly increase in response to favourable forcing that increases carbon uptake, the respiratory fluxes (R H and R A ) also increase and the rate of increase of F L declines.In the CLIM+CO 2 simulation, the per unit area sink of around 44 g C m −2 yr −1 (or equivalently 63 g C m −2 yr −1 over the forested area as shown in Fig. 5b) translates to an amount of ∼ 41 Tg C yr −1 when multiplied by the 944 700 km 2 area of the province.In Table 2 the magnitude of our simulated sink is compared to two observation-based estimates.The inversion-based estimate of net atmosphere-land CO 2 flux from Deng et al. (2007) for 2003 is 38 ± 66 g C m −2 yr −1 and the modelled value for the same year is 37 g C m −2 yr −1 .The caveats with modelled results are that about 6 % of the province's forest area was affected by the mountain pine beetle outbreak around that time and we do not account for losses due to harvest and insects, the resulting decay of dead organic matter or the rate of regrowth following those disturbances.In Table 2, our simulated sink of 63 g C m −2 yr −1 over the forested area in BC can also be compared to the NFCMARS-derived estimate for net ecosystem production (NEP) of 35 g C m −2 yr −1 (for the period 1990-2011 and updated from Stinson et al., 2011).However, the NFCMARS-derived estimates are based on the CBM-CFS3 model driven by empirical yield tables that do not take explicitly into account the effect of climate change and CO 2 fertilization, but the CBM-CFS3 does account for the impacts of all disturbances including harvest and insects.
The average sink enhancement over the forested area, in response to climate change and increasing CO 2 , of around 64 g C m −2 yr −1 for 1990-2005 may also be compared with the European forest carbon sink of 75 ± 20 g C m −2 yr −1 estimated by Luyssaert et al. (2010) over the same time period, given that both are located in the temperate mid-to highlatitude region.Of the cumulative realized sink of 2185 Pg C in the CLIM+CO 2 simulation, over the 1861-2000 period, about 63 % (1366 Pg C) is allocated to the vegetation pool ( H V ) and about 37 % (819 Pg C) to the litter and soil carbon pools ( H S ).Our distribution of the simulated sink into living (vegetation) and dead (litter and soil carbon) pools (63 % and 37 %, respectively) is comparable to estimates of 71 ± 15 % (woody biomass) and 29 ± 15 % (soil carbon) from Luyssaert et al. (2010) for European forests.
In Fig. 5 the simulated uptake is predominantly caused by changes in climate as opposed to CO 2 fertilization.Cumulatively over the 1861-2000 period, the sink caused by CO 2 fertilization is 519 Tg C and the sink caused by changing climate is 1623 Tg C, which combine almost linearly (2142 Tg C) to give about the same total carbon uptake as in the CLIM+CO 2 simulation (2185 Tg C).The weaker role of the CO 2 fertilization effect at high latitudes (the province of BC lies between 48 • and 60 • N) is consistent with the stronger CO 2 fertilization effect at low latitudes in the areas of higher temperatures (Hickler et al., 2008).Averaged over the period 1990-2005, the model simulates the BC-wide sink as 10.0 and 29.2 Tg C yr −1 in the CO 2 and CLIM simulations, respectively, and 42.3 Tg C yr −1 in the CLIM+CO 2 simulation (calculated using the provincial area of 944 700 km 2 ).
Although insect disturbances and harvest are not taken into account, these simulations do take into account the effect of fire.However, the area burned in BC is small.Averaged over the period 1970-2010 the annual forest area burned in BC is 726 km 2 (http://nfdp.ccfm.org/data/compendium/html/comp_31e.html),which translates into 0.08 % area burned annually.Over the same duration CTEM simulates a higher annual burn fraction of 0.28 %.The simulated area burned is modelled following the fire parameterization in CTEM (Arora and Boer, 2005b).The simulated annual burnt area and emissions (F F =∼ 0.8 Tg C yr −1 ) do not show any significant trend over the 1900-2009 period and therefore do not affect the simulated sink.Figure 5b shows forest-only net atmosphere-land CO 2 flux for forest regions in BC from the CLIM+CO 2 simulation.The averaged forest net atmosphere-land CO 2 flux over the 1990-2010 period is 63 g C m −2 yr −1 and averaged forest net primary productivity (NPP) is 355 g C m −2 yr −1 (Fig. 7b), which is lower than the NFCMARS-derived estimate of 447 g C m −2 yr −1 (updated after Stinson et al., 2011Stinson et al., , for the period 1990Stinson et al., to 2011)).The total sink enhancement from the CLIM+CO 2 simulation for forests only is ∼ 42 Tg C yr −1 (63 g C m −2 yr −1 multiplied by the area of the provincial forests, 663 981 km 2 ) The spatial distribution of model-simulated sink in response to changing climate and increasing CO 2 is shown in Fig. 6.The values peak over the centre of the province and over Vancouver Island.As expected, the sink occurs over areas with vegetation (see spatial distribution of LAI and GPP in Fig. 3).

Stem wood growth rate
Is the simulated response to changing climate and increasing atmospheric CO 2 concentration in our study realistic?Here we attempt to answer this by comparing our simulated stem wood growth to coastal BC ground-based measurements from Hember et al. (2012) in Fig. 7. Hember et  7a).Their data supported their hypothesis that the overall positive trends in stem wood growth rate were likely driven by long-term enhancement in NPP due to climate and atmospheric CO 2 forcing.We calculate stem wood growth averaged over coastal BC (essentially the grid cells with the coastal needleleaf evergreen PFT) using NPP allocated to the stem component in CTEM.Despite the inexact match in spatial coverage between estimates, CTEM-simulated (black and blue lines in Fig. 7a) multi-annual variability in stem growth is consistent with the observation-based estimate (green line in Fig. 7a).The inventory-based data reflect the aggregation of periodic (5 or 10 yr interval) measurements and therefore only represent variability on timescales exceeding the average interval length so the comparison is based on the 6 yr moving average CTEM values.The comparison of observation-based and simulated stem wood growth rate is therefore not completely consistent but nevertheless provides a valuable means of validating model simulations at this spatial and temporal scale.Observation-based data show an average stem wood growth rate of 3.45 Mg C ha −1 yr −1 over the period 1959-1998 and an average rate of increase of 1.5 % per year (green line in Fig. 7a).The simulated stem wood growth averages 2.8 Mg C ha −1 yr −1 and increases at a rate of 2.7 % per year (blue line in Fig. 7a).One possible reason for the higher observation-based average value is that the inventory data reflect productive forest sites, while no such distinction is made in CTEM, which represents an average over the forested fraction in a grid cell.It is also possible that the allocation to stem wood in CTEM is low.Higher values of NPP allocation to stem in the model, however, would imply even a stronger apparent increase in stem growth rates than we currently simulate.Both observation-based and simulated data show a decline in stem wood growth rate during the mid-1960s.During the 1980s, while the observation-based data show a decline, the simulated values show no trend.After 1990, both observation-based and simulated values show an increasing trend in stem wood growth rate.The simulated increase in stem wood growth in coastal and interior (not shown) regions is associated with an increasing trend in simulated NPP, which is shown in Fig. 7b.Despite inconsistencies in the observation-based and simulated stem wood growth, the comparison does highlight that the model captures the broad-scale increase in stem wood growth rate and the decadal trends over the 1959-1998 period that were observed in the forest inventory plots from coastal BC.

Discussion
There are at least two caveats associated with our results.The first is that CTEM does not take into account age-class distribution of trees, and the modelled response to changes in climate and CO 2 concentration is that of a tree of an average age in the landscape.In the real world, of course, ecosystems have irregular histories of disturbance and at the landscape scale there are trees of different ages.Hember et al. (2012) attempt to distinguish the effects of intrinsic (species composition, soil fertility and age class) and extrinsic (CO 2 fertilization and climate change) factors on stem wood growth.They use regression models to separate out the intrinsic factors (their Fig. 7b).In the province of BC, the forest inventory shows that average forest age is increasing and the agedependent reduction in tree growth counteracts the environmentally driven growth enhancement.Indeed when the effect of intrinsic factors is removed, the resulting rate of increase of stem wood growth attributed to extrinsic factors is 3.0 % (red line in Fig. 7a), bringing it into closer agreement with the CTEM-simulated rate of increase of stem wood growth of 2.7 %.
The second caveat is that CTEM does not model harvesting of wood and insect disturbances.Our BC-wide estimate of sink of 44 g C m −2 yr −1 over the 1980s, 1990s and 2000s (or equivalently 63 g C m −2 yr −1 over the forested area) is the response of BC's terrestrial ecosystems to changing climate and increasing CO 2 .However, from a greenhouse gas accounting perspective, where the harvested wood is removed from the forests, the carbon source generated by harvesting, the decay of the resulting slash and the subsequent forest regrowth also needs to be taken into account.The carbon source from harvesting may not be generated in the province of BC itself since most forest products are exported.Stinson et al. (2011), who estimate the carbon budget over Canada's forests using the CBM-CFS3, provide an estimate of 30 g C m −2 yr −1 for harvest-related losses.This estimate, however, does not include the effects of decaying post-harvest residues and regrowth of the harvested forest although these processes are included in the CBM-CFS3.
Since 2001 the province of BC has been affected by the mountain pine beetle insect disturbance, which has killed a large fraction of forests in the interior of the province (Kurz et al., 2008;Stinson et al., 2011), yielding a large carbon source to the atmosphere.Our historical simulations do not include this disturbance and the simulated sink after 2001 is thus overestimated.Our simulations also do not include other prior insect disturbances in the province (e.g.see Kurz et al., 1996).Insect disturbances cause foliage loss and/or tree mortality and lead to reduced carbon uptake and increased loss due to tree mortality from forest ecosystem during outbreak periods (Kurz et al., 2008).For the peak years of the mountain pine beetle outbreak in BC, the reduction in carbon uptake due to tree mortality and the increase in R H from beetlekilled trees lowered the annual carbon sink by 20 Tg C yr −1 (Kurz et al., 2008, their Fig. 4), which translates to an additional source of about 30 g C m −2 yr −1 when divided by the 663 981 km 2 forest area in our simulation.
Carbon budget assessments for forests routinely take into account harvesting, insect and fire disturbance losses.These assessments are, however, based on models that simulate tree growth using empirical yield curves; for example, the CBM-CFS family of models (Kurz et al., 2009) that is used by the Canadian Forest Service (Stinson et al., 2011, and references therein).Tree growth in the CBM-CFS models does not explicitly account for changes in tree growth due to changes in climate and CO 2 fertilization.The effect of these processes is included in so far as these processes have affected past tree growth and their effect is reflected in empirical yield curves.Our results suggest that changing climate and increasing CO 2 have both contributed to an increased carbon sink for the province of BC, and inclusion of these processes in forest carbon budget assessments will lead to higher estimates of carbon uptake.
The results reported here are dependent on the climate data used to drive the models used in this study.There are limitations in the driving CRUNCEP climate data.In particular it does not show a large positive precipitation trend over BC, like the station-based data.Photosynthesis at high latitudes is generally not soil moisture limited except during the peak summer months.Had the CRUNCEP data exhibited a large positive trend in precipitation, like the station-based data, we would expect the simulated sink to be somewhat stronger.

Summary and conclusions
Several lines of observation-and modelling-based studies suggest that the mid-to high-latitude ecosystems are currently a sink of carbon.Here, we have used a processbased ecosystem model to simulate the response of terrestrial ecosystems in the province of BC, Canada, to changing climate and increasing CO 2 over the historical 1860-2010 period.Our results suggest that these two forcings have resulted in an enhancement of BC's terrestrial ecosystems sink of carbon, relative to pre-industrial conditions, of around 44 g C m −2 yr −1 (or 63 g C m −2 yr −1 over BC's forests) during the 1980s, 1990s and 2000s.
To assess if the simulated response to changing climate and increasing CO 2 in our study is realistic we have compared the simulated stem wood growth to observationbased inventory estimates for coastal BC from Hember et al. (2012).Our rate of increase of stem wood growth of 2.7 % over the 1959-1998 period compares well with the observation-based estimate of 3.0 % from Hember et al. (2012) when the effect of intrinsic factors (including age class, soil fertility and species composition) is accounted for.This comparison corroborates the combined model response to changing climate and increasing CO 2 .The model validation against stem wood growth rates from inventory plots from coastal BC provides confidence in the simulated response of the province's ecosystems to changing climate and increasing atmospheric CO 2 concentration.This validation also provides the groundwork required to use our modelling framework for investigating the response of the province's terrestrial ecosystems to future changes in climate and atmospheric CO 2 concentration.
The per unit area sink enhancement estimate of 44 g C m −2 yr −1 during the 1980s, 1990s and early 2000s translates to an amount of 41 Tg C yr −1 .Including harvestrelated losses including those from burning or decay of postharvest residues and from the recent mountain pine beetle infestation will further reduce the sink estimate.The estimated sink enhancement suggests that current estimates of the net forest carbon balance of the province, based on empirical yield curves, may be underestimated because they do not fully account for the growth enhancements due to changes in climate and increase in atmospheric CO 2 concentration.
Although limited in their domain, our regional-scale results contribute to the existing evidence that mid-to highlatitude terrestrial ecosystems are currently responding to changing climate and increasing CO 2 by acting as a sink www.biogeosciences.net/11/635/2014/Biogeosciences, 11, 635-649, 2014 of atmospheric carbon.Specifically, our study estimates that about three-quarters of the simulated historic sink enhancement in the province of BC is attributable to changing climate.

Appendix A
Dynamic global vegetation models (DGVMs) typically use vegetation classification based on the PFT concept, and species-level differences are not considered.The PFT concept has emerged from the school of thought that suggests vegetation classification with a combined ecological and biophysical focus is possible in relation to a plant's form and the manner in which it interacts with its environment (Box, 1996).CTEM's use of a single needleleaf evergreen PFT works reasonably well at the global scale (e.g.Arora et al., 2009) in an Earth system model framework where grid resolutions are typically 200-300 km.However, for application over BC at a spatial resolution of ∼ 40 km, we found that, while this single needleleaf evergreen PFT yields reasonable LAI and GPP in the coastal region compared to observationbased estimates, it yields unrealistically low LAI and vegetation biomass in the interior of BC (left panels in Fig. A1 compared to right panels in Fig. 3).As a result, the simulated province-wide averaged GPP of 398 g C m −2 yr −1 for the period 1998-2005 is low compared to the observation-based estimate of 597 g C m −2 yr −1 based on Beer et al. (2010) shown in Fig. 3.In Table A1, for the single needleleaf evergreen PFT case, the province-wide areally averaged summer maximum LAI of 1.4 m 2 m −2 is also low compared to observation-based estimate of 2.5 m 2 m −2 from Deng et al. (2006) (see also Table 2 and Fig. 3).
The BC Ministry of Forests, Lands and Natural Resource Operations (BC MFLNRO) divides the province into 16 biogeoclimatic zones characterized mainly by the climate and the species of dominating conifers.Coastal BC is occupied primarily by western hemlock, western red cedar and coastal Douglas fir, while other needleleaf evergreen species (primarily pines, spruces, and subalpine fir, but also interior Douglas fir, western hemlock and western red cedar) occupy the interior of the province (BC Ministry of Forests and Range, 2008).Being a DGVM, CTEM does not represent the level of detail necessary to parameterize species-level differences, nor are there data easily available which may be used to estimate CTEM parameter values that would reflect specieslevel differences.For this modelling study, therefore, we make the first-order distinction based on climate and split CTEM's default needleleaf evergreen PFT into coastal and interior types.The coastal type retains the default parameter values of CTEM's needleleaf evergreen PFT, including a leaf life span of 1.5 yr.For the interior needleleaf evergreen PFT, we assume 50 % lower rates of leaf loss from cold and drought in the phenology parameterzation of CTEM (Arora and Boer, 2005a), consistent with the colder and drier cli- mate in the interior of the province to which the trees have adapted.In addition, we assume a longer leaf life span of 5 yr for the interior type following Reich et al. (1995), who suggest that greater conifer leaf longevity in colder environments is associated with thicker leaf cuticles and other structural modifications to minimize winter desiccation.Following the BC MFLNRO's biogeoclimatic zones, we assign the coniferous forests west of the coastal mountains as CTEM's coastal needleleaf evergreen PFT and those to the east of the coastal mountains as interior needleleaf evergreen PFT.
In Fig. A1, simulated LAI and GPP increase in the interior of the province for the two needleleaf evergreen PFT simulations compared to the simulation with CTEM's default needleleaf evergreen PFT yielding better agreement with observation-based estimates (Fig. 3).Table A1 shows that province-wide averaged GPP and LAI, and provincewide total terrestrial pools of vegetation, litter and soil carbon, are higher for the simulation with two needleleaf evergreen PFTs and compare better with available observationbased estimates in Table 2.

Copyright statement
The works published in this journal are distributed under the Creative Commons Attribution 3.0 License.This license does not affect the Crown copyright work, which is re-usable under the Open Government Licence (OGL).The Creative Commons Attribution 3.0 License and the OGL are interoperable and do not conflict with, reduce or limit each other.© Crown copyright 2014

Fig. 1 .
Fig. 1.Spatial distribution of fractional vegetation cover for all PFTs and for needleleaf evergreen trees, broadleaf cold deciduous trees and C 3 grass PFTs.
shows the trends in BC province-wide averaged climate variables in the CRUNCEP data that are used to drive CLASS and CTEM.The values for pressure, wind speed and longwave radiation are prescribed constant from 1901 to around 1940 in the CRUNCEP data over the BC region likely because of unavailability of observation-based estimates.Over the 1901-2010 period, all province-wide averaged climate variables show positive trends except radiation.The province-wide averaged temperature has increased by about 1.0 • C. The province-wide averaged precipitation shows a weak positive trend of 4.4 % increase over the 1901-2010 period.The CRUNCEP data, however, have limitations.
averaged over the 2000-2005 period.The observation-based estimate of GPP is from Beer et al. (2010), who analyse the ground-based carbon flux tower observations from about 250 stations and apply diagnostic models to extrapolate them to the global scale.This GPP estimate corresponds to the period 1998-2005.Being a global product, the GPP estimates from Beer et al. (2010) cannot represent the detailed spatial patterns for the province of BC but nevertheless provide an observation-based estimate.Three different observation-based estimates of vegetation biomass are used.The first observation-based vegetation biomass estimate for BC is based on a global data set of vegetation carbon stocks in the year 2000 and uses a range of spatially explicit climate and vegetation data sets to map vegetation biomass values

Fig. 2 .
Fig. 2. BC province-wide averaged climate variables from the CRUNCEP data that are used to drive CLASS and CTEM models for the period 1901 to 2010.The dark-red curves are 10 yr moving averages and the blue lines are linear trends whose values are also noted.The trend values are normalized to change per 100 yr.

Fig. 3 .
Fig. 3. Geographic distribution of simulated maximum leaf area index (LAI) and gross primary productivity (GPP) for BC compared to observation-based estimates.The left column shows results from the CLIM+CO 2 simulation that considers coastal and interior needleleaf evergreen PFTs separately.LAI observations are from Deng et al. (2006) averaged over the 2000-2005 period, and GPP observations are from the global product of Beer et al. (2010).The model results are averaged over the same respective periods as well.

Fig. 4 .
Fig. 4. Spatial distribution of simulated vegetation biomass (left panel) and litter and soil carbon mass (right panel) (Kg C m −2 ).Data are derived from the CLIM+CO 2 simulation and averaged over the period 1990-2005.

Fig. 5 .
Fig. 5. (a) BC-wide simulated net atmosphere-land CO 2 flux (gC m −2 yr −1 ), in response to changing climate and/or CO 2 concentration, during the 1900-2010 period for three historical simulations -CO 2 , CLIM and CLIM+CO 2 -as explained in Sect.2.4.The lines are 10 yr moving averages.Positive values indicate a carbon sink over land, and negative values imply a source of carbon to the atmosphere.(b) Forest-only net atmosphere-land CO 2 flux compared to province-wide averaged net atmosphere-land CO 2 flux from the CLIM+CO 2 simulation, excluding the impacts of forest insect disturbances and harvesting.

Fig. 6 .
Fig. 6.Annual mean simulated net atmosphere-land CO 2 flux, in response to changing climate and CO 2 concentration for the period 1990-2005.Results are from the CLIM+CO 2 simulation and units are gC m −2 yr −1 .Positive values indicate a carbon sink over land, and negative values imply a source of carbon to the atmosphere.

Fig. 7 .
Fig. 7. Comparison of simulated and ground-based observed stem wood growth rate of conifer trees in coastal BC.(a) Observations are based on data collected over the period 1950-2002 as explained in Hember et al. (2012).Model estimates are for all needleleaf evergreen trees exclusively in coastal region of BC.Black and blue lines for simulated results are 6 yr moving averages.Linear regression lines are on top of the model estimates and observations during the 40 yr period for which observations are available.Panel (b) shows forest-only net primary productivity (NPP) in all forest-covered regions in BC, as well as 6 yr moving average of forest NPP.

Fig. A1 .
Fig. A1.Geographic distribution of simulated maximum leaf area index (LAI) and gross primary productivity (GPP) for BC.The left column shows results from the CLIM+CO 2 simulation with CTEM's default needleleaf evergreen PFT, and the right column from the CLIM+CO 2 simulation that considers coastal and interior needleleaf evergreen PFTs separately.The model results are averaged over the reference 1990-2005 period.

Table A1 .
Comparison of British Columbia total carbon pools and spatially averaged carbon fluxes averaged over the reference period 1990-2005 from the CLIM+CO 2 simulation for cases where one and two needleleaf evergreen PFTs are considered.Simulation with one needle Simulation with two needle leaf evergreen Net primary productivity (g C m −2 yr −1 ) 202 294 Gross primary productivity (g C m −2 yr −1 ) 419 633

Table 1 .
Plant functional types (PFTs) for which terrestrial ecosystem processes are modelled by CTEM in its default configuration at the global scale and for this study.CTEM PFTs that are present in the province of BC (see Fig.1) are shown in italic.

Table 2 .
Comparison of simulated vegetation biomass, leaf area index, gross primary productivity, net atmosphere-land CO 2 flux and stem wood growth rates with observation-based estimates for the province of British Columbia, Canada, from the CLIM+CO 2 simulation with the modified PFT set.Modelled results are reported for the 1990-2005 reference period but also for the time period that corresponds to each observation-based estimate.

Coastal region averaged stem growth rate
a Both the Penner et al. (1997) estimate and the updated estimate from Stinson et al. (2011) are for aboveground biomass.b Model estimates are for the 1990-2010 period because model simulations stop in 2010.