Study area
The pilot study area is located in the boreal forest of northern Alberta in the vicinity of the city of Ft. McMurray (Fig.1). The region is subject to a wide range of natural and anthropogenic disturbances, including intensive oil and gas exploration and development activities, forest harvesting, wildfire and insect outbreaks (Fig. 1, [40,41,42,43,44,45]). The pilot study area covers 2,482,770 ha where 1,267,725 ha are uplands (the focus of this study) with Luvisolic, Brunisolic and Gleysolic soils and 787,361 ha are wetlands with Organic peatland soils [46,47]. It is in the Central Mixedwood subregion of the Boreal Mixedwood forest of northern Alberta that has a subhumid to subarid (annual precipitation of 389 mm), cool continental climate (mean annual temperature 1.5° C), with long cold winters and warm summers [46]. In the Central Mixedwood subregion the dominant tree species is trembling aspen (Populus tremuloides Michx.) and co-dominant species are balsam poplar (Populus balsamifera L.), black spruce (Picea mariana (Mill.) BSP), white spruce (Picea glauca (Moench) Voss) and jack pine (Pinus banksiana Lamb.) [46].
The location of the pilot study area was chosen to represent the suite of disturbance types in the OSR, where the best data layers were available for annual change in disturbances over time, fine-scale disturbances from oil and gas exploration and development, and forest inventory with associated yield curves. The data layers were used to generate spatially-explicit inputs for the GCBM.
Modelling with the GCBM
Modelling of upland forest types was conducted using the GCBM, the next-generation, spatially-explicit version of the CBM-CFS3. The GCBM is composed of science modules linked to the Full Lands Integration Tool (FLINT) framework, an open source, modular, spatially-explicit modelling framework developed jointly by experts from Australia (Mullion Group), Canada (Canadian Forest Service) and the moja global organization [48] (Fig. 2). The science modules in the version of the GCBM used here replicate the representation of C science of the CBM-CFS3 but are linked to the FLINT, which assists in the processing of the spatial information. The CBM-CFS3 simulates the entire landscape one year at a time. In contrast, the GCBM simulates each pixel over the entire time series before moving on to the next pixel, which enables the use of parallel processing computing architecture and allows the spatially-explicit application of the GCBM to much larger landscapes. However, this approach requires that disturbance (and forest management) information is provided as input in spatial layers that define both the year and type of disturbances.
The GCBM uses, as inputs, a combination of a spatial forest inventory (including information like the location of forest types, tree species, stand age), mean annual temperature, and annual disturbance layers, along with non-spatial parameters or parameters that are spatially-referenced to terrestrial ecozones or administrative regions [29]. Model parameters include yield curves, volume-to-biomass conversion coefficients, and disturbance matrices to estimate the annual C balance of a study area. Spatial data for the GCBM are stored in a custom format using the EPSG 4326 / WGS 84 projection [49]. A Python tool is used to crop, re-project, and re-sample spatial layers from their original raster or vector format into the format and projection used by the model. Simulations in this project were conducted at 30 metre resolutions. Non-spatial data are stored in an SQLite database populated by a tool that loads a user-provided set of yield curves and a file defining any transitions of stand characteristics following disturbance, and imports other default model parameters from the CBM-CFS3 input databases. In this pilot study the only transitions following disturbance that were represented were reductions in stand growth following seismic line disturbances. Growth for the pixel after disturbance was reduced by the same proportions used for disturbed area for the old (0.27, i.e., established 1960 to 1999) and modern (0.17) seismic lines (see Results section “Energy sector disturbances” for rationale).
Disturbance matrices
The effect of different disturbance types on forest C in the year of the disturbance are modelled in the GCBM (and CBM-CFS3) using disturbance matrices (DMs) [28]. The DMs define the proportions of C that are transferred between model pools, between pools and the atmosphere, and between pools and harvested wood products, at the time of the disturbance. Default DMs used in the GCBM are well developed for wildfires, harvesting, and insect outbreaks [29, 30]. Currently there are two default DMs in the GCBM that could be used to represent disturbances from oil and gas exploration and development. They are intended for application to situations that invoke a land-use change that typically would be associated with large open-pit mining operations used for extraction of resources such as gold, copper, aluminum and coal. The default DMs could have been applied to surface mining in this pilot study, but new DMs were developed in this project for all oil and gas disturbance types, including surface mining, to take advantage of the most recent understanding of the impacts of these disturbance types specific to the OSR.
A team of scientists with expertise in reclamation field research in the OSR and those with expertise in forest C accounting and reporting, and boreal forest C modelling, participated in a two-day workshop to define the types of DMs required for modelling the effects of oil and gas exploration and development on forest C in the oil sands region of Alberta. Scientists attending are leaders in the fields of surface mining reclamation, reclamation land disturbed by in-situ development and exploration disturbances. Others were experts in modeling using the CBM and GCBM and its appropriate implementation for reporting purposes. Initially, discussions were held to identify the suite of oil and gas related disturbance types in the region and factors that would affect fate of woody materials (i.e., proximity to a mill). These disturbance types were distilled into a list for which expert opinion was that sufficient experience and understanding existed to estimate the proportions for C transfers in a disturbance matrix. These new DMs were then populated using consensus-based values quantifying the proportion of C transfers. Each expert independently estimated the proportion for each transfer in each disturbance matrix. Values were accepted where the same proportion was estimated by all experts. A discussion was held for each proportion where experts’ initial estimates differed, in order to arrive at a consensus. The intent was to create a suite of DMs for disturbances types that can be identified using the available spatial data layers, and for which there is sufficient understanding of C transfers in response to the disturbances. The DMs were developed for upland situations only, and took into consideration factors such as ecosite type, distance from a mill, permanence, and vintage of seismic lines (see Results section “Energy sector disturbances” for details).
Model inputs
Ecological parameters for annual processes
Annual processes (e.g., decomposition, mortality) determine transfers of C from living biomass to standing and downed deadwood pools, from biomass and deadwood to organic and mineral soil pools, and to the atmosphere are explicitly simulated in the GCBM in the same way as in the CBM-CFS3 (Fig. 3, [28]). Carbon is physically transferred from live to dead pools because of annual biomass turnover and stand mortality (i.e., a yield curve with declining volume). Disturbances transfer C from biomass to dead organic matter (DOM) pools, to the atmosphere, and to harvested wood products. Carbon is also transferred from the slowest decay pool in the organic soil horizons to the slowest decaying pool in the mineral soil, representing a transfer as dissolved organic C. Carbon is moved between DOM pools and to the atmosphere as a result of decomposition. Base decay rates in the GCBM are modified in response to mean annual temperature using a Q10 relationship. Base decay rates and Q10 vary by DOM pool. Parameters used in this study are those specific to the Boreal Plains ecozone [51] and are specified in Kull et al. [29]. In this study we assume that all C transferred to the forest product sector is instantly oxidized and released to the atmosphere, and thus overestimate the direct emissions from harvest by not estimating C retention in harvested wood products and landfills [30, 31, 34].
Inventory data layers
This project integrated an Alberta Vegetation Inventory (AVI, [52]) forest inventory provided by Alberta Pacific Forest Industries, Inc. (AlPac) with two strata defining the softwood and hardwood component. The stratum raster was created by rasterizing the AVI to 30m spatial resolution and filling the null pixels with species from a basemap derived by the Canadian Centre for Remote Sensing (CCRS) land cover time-series [44]. The basemap was created using the values of the pixels from the time-series that had a forest cover for the first three years of the time-series. The conifer, deciduous and mixedwood forest classes were then converted to the stratum type that corresponded to the dominant AVI stratum intersecting the forest class. To match the yield curve stratification, conifers were labelled PjO-CFM (describing open and closed jack pine led conifer stands located on fair and medium productivity sites), mixedwoods labelled as AwS_N (referring to aspen (Populus spp.) and spruce (Picea spp.) mixedwoods stands located on the northern portion of the AlPac's forest management agreement area), and deciduous labelled Aw_comp (comprising of an all pure aspen stand).
Using AlPac’s AVI attributes, every mapped polygon in the AlPac landbase within the pilot study area was post-stratified into one of the 22 yield strata that separate stands or stand groups with different growth characteristics. The key characteristics used to stratify the landbase included species composition, crown closure, timber productivity rating, understory occurrence, and geographical location [53]. The yield strata assignment was computed using a SAS script incorporating rules used to define yield strata (rules taken from Table 10 in [53]).
Yield curves
The yield strata were used to link pilot study inventory with yield curves. Yield strata assignments involved the determination of stand characteristics such as overstory and understory broad cover types or leading conifer species [53]. AlPac’s detailed Forest Management Plan is supported by a timber supply analysis, which entails development of yield curves for the various AVI strata. For each yield stratum, defined by AVI stand attributes including species composition, crown closure, site productivity and location (Tables 10 and 12 in [53]), a series of empirical yield curves was constructed based on the AlPac’s temporary sample plot and permanent sample plot data. Merchantable yield curves were fit using non-linear regression, applied separately to the total, coniferous, and deciduous, volume-age pairs [53]. Coniferous and deciduous curves generated separately in each yield stratum were matched to species in GCBM. Merchantable yield estimates were projected using five-year intervals, which were converted to annual estimates using linear interpolation.
Disturbance data layers
A merged disturbances time-series (1985 - 2012) was created from the CCRS disturbance time-series [44] and the 2012 Alberta Biodiversity Monitoring Institute (ABMI) human footprint [54]. Year and type of disturbance for wildfire, insect disturbances, harvest, and surface mines were taken directly from the CCRS time-series. For other disturbance types, the ABMI time-series was used to assign disturbance type, and the CCRS time series was used to assign a year to that disturbance type because, at the time this work was completed, ABMI did not provide disturbance year attributes. The CCRS identified approximately 1,602 ha of generic disturbances with/without regeneration that included the year of disturbance. These were generally related to oil and gas activities. To assign a year to the ABMI disturbance events, for each ABMI disturbance object, we extracted the mode from the year of first disturbance of the object in the CCRS time series. For well pads, pipelines, roads, residential and industrial sites in ABMI, only the disturbances intersecting with at least one CCRS pixel were included in the merged disturbance times-series. For seismic lines, we started by adding the disturbance events that intersected at least one CCRS pixel. We then used the disturbance year of the nearest neighbor to add year attributions to the remaining seismic lines that did not intersect with CCRS pixels. All the CCRS generic disturbances pixels that did not intersect with ABMI disturbances were retained. Disturbances prior to 1985 were excluded from the time series, so effects from surface mining prior to 1985 were not quantified. Seismic line disturbances are narrower than the width of a pixel (approximately 30 m) and were represented by reducing the total disturbance effect of a seismic line type proportional to the disturbed area of the pixel (see RESULTS section “Energy sector disturbances”).