Stand-clearing disturbances, which remove most of the tree cover but are followed by forest regrowth, affect extensive areas annually, yet each event is usually much smaller than a typical grid cell in Earth system climate models. This study argues that the approach taken to account for the resulting subgrid cell dynamic heterogeneity substantially affects the computation of land–atmosphere exchanges. The authors investigated in a simplified model the effects of three such approaches on the computation of albedo over boreal forests. It was found that the simplest approach—in which any new disturbance-created patch was immediately merged with the rest of the grid cell—underestimated the annual reflected solar radiation by ~3 W m−2 on average (a relative error of 15%) compared with the most accurate approach—in which albedo computations were performed for each individual subgrid patch. This study also investigated an intermediate approach, in which each patch was tracked individually, but albedo was estimated from a much smaller number of subgrid tiles grouping patches having a similar amount of tree cover. Results from this third approach converged quickly toward the most accurate results as the number of tiles increased and were robust to changes in the thresholds used to assign patches to specific tiles. When computing time prevents implementing the most accurate approach in Earth system climate models, the results advocate for using strategies similar to the intermediate approach in order to avoid biasing the net radiative forcing of stand-clearing disturbances toward a warming impact, at least over boreal forests.
Each grid cell in coupled climate–terrestrial vegetation models is bound to include major spatial heterogeneities that span various orders of magnitude and influence the processes represented (Giorgi and Avissar 1997). The basic approaches accounting for subgrid cell heterogeneity in these models can be divided into two main categories. The composite (also named aggregated) approach computes land–atmosphere exchanges as a function of a single “representative” state of the grid cell (Giorgi and Avissar 1997). The approach allows for the coexistence of different plant functional types (PFT), but these all experience the same environmental conditions. The second category groups the mosaic, mixture, and statistical–dynamical approaches, which all account explicitly for the contribution of different portions of the grid cell. The mosaic approach subdivides the grid cell into tiles that respond independently to the same climate input variables so that soil temperature, for example, can differ among tiles (Avissar and Pielke 1989). Different criteria can serve to subdivide the grid cell into different tiles, and the approach is often applied in a nested framework: for example, first to separate vegetated from nonvegetated areas (e.g., lakes) and then to subdivide the vegetated area into different PFT tiles (Bonan et al. 2002). Given that each tile is homogeneous, the mosaic approach implicitly uses the composite approach, albeit at a finer scale (i.e., within a PFT tile instead of across the entire grid cell). The mixture approach is a variation of the mosaic approach, in which the latent and sensible heat fluxes from the different tiles of the same grid cell interact together (Koster and Suarez 1992). The statistical–dynamical approach aims to represent the full range of heterogeneity across the entire grid cell through continuous probability density functions (Avissar 1992). The resulting equations can be integrated analytically for very simplified representations of land–atmosphere interactions, but applying the approach to state-of-the-art Earth system climate models seems hardly feasible because of the practical (discretization of the probability density functions) and conceptual (modification of each probability density function in response to climate change) issues. In large-scale climate models, which currently do not resolve mesoscale subgrid circulation (Pielke et al. 2007), none of the previous approaches account for the specific location of the heterogeneous elements within the grid cell.
The previous approaches have been developed to represent static subgrid cell heterogeneity (Giorgi and Avissar 1997). Research on how to better account for the impacts caused by static heterogeneity is still ongoing, particularly in the context of topographic variations (Ke et al. 2013; Newman et al. 2014; Rotach et al. 2014; Subin et al. 2014). Yet the spatial heterogeneity of the land surface is dynamic, among others due to changes in vegetation cover resulting from human or natural causes. Over the last centuries, agriculture has led to substantial and permanent (or long term) changes in land cover (Ramankutty and Foley 1999). These deforestation events can adequately be simulated in climate models by permanently decreasing the grid cell fraction occupied by the forested tile(s) and increasing the agricultural tile(s) (e.g., Matthews et al. 2004). (Afforestation is more difficult to represent because the regrowing forest differs considerably from the rest of the forest for many years.) Other anthropogenic and natural disturbances, however, lead to more complex changes in spatial heterogeneity through the temporary removal of tree cover over extensive areas, followed by forest regrowth. Fire and logging are important examples of such events, which we define as stand-clearing disturbances. At the global scale, fire burns about 350 Mha yr−1 (Giglio et al. 2013) and emits 1.5–3 PgC yr−1 (Mieville et al. 2010; van der Werf et al. 2010), while logging impacts about 30 Mha yr−1 and removes between 0.2 and 0.9 PgC yr−1 from forests (Hurtt et al. 2006; Pan et al. 2011). Although the previous numbers include some treeless ecosystems and permanent changes (e.g., deforestation fires), fire and logging events that are followed by tree regrowth affect extensive areas of forest each year. These events have been shown to be consequential for climate through their impacts on carbon cycling and energy exchanges, including the possibly large effect of aerosols in the case of fire (Jacobson 2004; Jones et al. 2007; Unger et al. 2010; Lawrence et al. 2012; Ward et al. 2012; Jacobson 2014; Landry et al. 2015). Furthermore, stand-clearing disturbances play a major role in ecological succession and can trigger vegetation shifts (Turner 2010). For all these reasons, representing the transient impacts from stand-clearing disturbances would improve the simulation of vegetation dynamics, land surface properties, and exchanges of carbon, energy, and water with the atmosphere.
Stand-clearing disturbances introduce dynamic (i.e., year after year) changes in subgrid cell heterogeneity, even if the area of the forested portion of the grid cell stays constant. Representing these recurring events is thus challenging, particularly for the first generation of Dynamic Global Vegetation Models (DGVM). These DGVM resort to large-scale parameterizations that, although usually distinguishing various PFT, do not explicitly treat subgrid heterogeneity for the same PFT (Prentice et al. 2007; Quillet et al. 2010). Examples of first-generation DGVM include IBIS (Foley et al. 1996; Kucharik et al. 2000), TRIFFID (Cox 2001), LPJ-DGVM (Sitch et al. 2003), CTEM (Arora and Boer 2005), ORCHIDEE (Krinner et al. 2005), JSBACH (Raddatz et al. 2007; Brovkin et al. 2009), and CLM4 Oleson et al. 2010) (expansions of acronyms are available at http://www.ametsoc.org/PubsAcronymList). Second-generation DGVM are more amenable to the representation of subgrid heterogeneity because they simulate small-scale interactions among many individual plants (including for the same PFT) during the establishment, growth, and mortality phases (Fisher et al. 2010). Most of these second-generation DGVM consist of stochastic gap models and include Hybrid (Friend et al. 1997), LPJ-GUESS (Smith et al. 2001), and SEIB-DGVM (Sato et al. 2007). The other class of second-generation DGVM are composed of the ecosystem demography (ED) family, which simulates the mean behavior of an ensemble of stochastic individual-based models over a large spatial scale through a deterministic “scaling up” approximation (Moorcroft et al. 2001; Medvigy et al. 2009). Second-generation DGVM are, however, considerably less computationally efficient than large-scale parameterizations (Scherstjanoi et al. 2013), which become problematic for large-scale and long-term climate simulations. Consequently, the vast majority of Earth system climate models currently use first-generation DGVM.
Various efforts have recently been made to “scale down” first-generation DGVM, usually aiming to better represent demographic processes as well as the subgrid vegetation heterogeneity caused by land-use changes (Shevliakova et al. 2009), forest self-thinning, human-caused thinning, and stand-clearing harvest (Bellassen et al. 2010), generic stand-clearing disturbances (Scherstjanoi et al. 2013), or fire (Haverd et al. 2013; Yue et al. 2013; Haverd et al. 2014). These studies focused on improving the simulation of the carbon cycle directly, or indirectly through more realistic forest structural attributes (e.g., allometry), but did not consider the biogeophysical fluxes that also influence the climate. Reduced evapotranspiration and surface roughness from deforestation increases the local temperature almost everywhere on land, but the cooling impact from the increased albedo has a stronger effect globally and over most extratropical regions (Davin and de Noblet-Ducoudré 2010; Lee et al. 2011). In boreal regions, the albedo-induced cooling (warming) caused by permanent forest cover reduction (increase) has generally been found to be greater than the associated warming (cooling) from reduced (increased) terrestrial carbon storage (Betts 2000; Claussen et al. 2001; Bala et al. 2007; Bathiany et al. 2010; Bernier et al. 2011) or to at least have the same magnitude (Arora and Montenegro 2011; Pongratz et al. 2011). In the case of stand-clearing fire in boreal forests, studies have also concluded that the albedo cooling was stronger than the carbon warming (Randerson et al. 2006; O’Halloran et al. 2012).
The purpose of our study was to quantify how different representations of the subgrid cell dynamic heterogeneity caused by fire or logging stand-clearing disturbances affect the modeling of albedo in boreal forests. In section 2, we present the Heterogeneous Landscape Model (HLM) we developed for this purpose as well as the three generic approaches we assessed. Using the simpler HLM instead of state-of-the-art models facilitated the interpretation of results (section 3) and allowed us to obtain more generalizable outcomes. In section 4, we discuss the importance of accounting for the impacts of subgrid dynamic heterogeneity in climate models as well as the implications of the HLM results on the representation of this heterogeneity in existing models.
2.1. Heterogeneous Landscape Model
The HLM reproduced the computation of solar radiation exchanges used by many climate models (i.e., as the mean value of fluxes from the vegetated and “canopy gap” portions of the grid cell), but the HLM used empirical data for the different variables required instead of estimating them through prognostic simulations. We modeled a grid cell representative of the black spruce [Picea mariana (Mill.) B.S.P.] boreal forest in Canada, with albedo values that were the same from one year to another but did vary daily. We compared the grid-level mean annual reflected solar radiation (, in W m−2) estimated by three approaches accounting differently for the effect of subgrid cell dynamic heterogeneity. We assumed that the total vegetation carbon in the grid cell was the same for the three approaches (an assumption easily implemented in the HLM due to its simplified dynamics, but that would have been difficult to impose in state-of-the-art models), so that differences in RSOL among approaches came only from the manner in which they accounted for the subgrid distribution of this carbon.
Spatial heterogeneity in the grid cell came from stand-clearing disturbances only, which were either logging or fire. Each year, a stand-clearing disturbance event affected a constant proportion of the grid cell (unitless) and removed all trees from the patch thus created, after which even-aged trees grew back in the patch. We assumed that the grid cell was in equilibrium with a specific disturbance regime; hence, the age-class distribution of patches remained constant, although the state of each individual patch kept changing from one year to another. Logging and fire varied in the spatial distribution of disturbance events: logging happened in the oldest patch only, whereas fire affected all patches in proportion to their area. Consequently, logging resulted in a number of same-sized patches that was equal to the rotation length (i.e., the number of years before a patch is logged again), with the youngest patch aged zero and the oldest patch aged one year less than the rotation length. Fire, on the other hand, led to an area-based age-class distribution given by ()t (Van Wagner 1978), where t is the patch age (an integer varying formally from 0 to , in years). We made the previous choices for modeling purposes, not to accurately represent actual disturbance regimes.
2.1.2. Main equations
RSOL was the mean value of the daily amount of reflected solar radiation, over the 365 days (i) of the year:
where sol is the daily incoming solar radiation at the ground level (W m−2), () is the mean daily albedo for tree-covered (treeless) areas (unitless), and is the fraction of the total grid cell covered by trees (unitless). Only varied among the three approaches assessed, but it was always based on the function estimating the tree cover fraction (unitless) over a given area. Existing models compute in different manners, but generally based on some carbon-related variables (e.g., leaf area index). We used the function from IBIS (Kucharik et al. 2000), which computes directly from the amount of stem (including branches) carbon (, in kgC m−2):
Together, Equations (1)–(2) constituted a generic computation of grid-mean surface albedo (or RSOL) as the weighted average of tree-covered and canopy gap , with the weight () being computed as a function of in a manner that depended upon the specific approach (see section 2.1.3 below). These two equations therefore involved a nonlinear dependence of RSOL on —which is realistic because vegetation cover saturates with increasing biomass—that ultimately explained why albedo differed among the three representations of subgrid cell dynamic heterogeneity. The yearly dynamics of was given by
where a is the proportion (unitless) of net primary productivity (NPP, in kg C m−2 yr−1) going to stem, and f is a fraction (unitless) that is equal to 1.0 minus turnover rate. Note that the first term on the right-hand side of Equation (3) is implicitly multiplied by yr and that . The different values for , , , , a, and f required by the HLM (see section 2.1.4 below) came from empirical data, mostly representative of black spruce boreal forest in Canada.
2.1.3. Approaches assessed
The three approaches could all, conceptually at least, be implemented into real climate models, so we first present them in general terms before providing more details regarding their implementation in the HLM. The first approach consisted of adding one specific tile for each new patch created by a stand-clearing disturbance event and keeping track of all patches indefinitely. This Full approach would provide the most accurate results regarding subgrid dynamic heterogeneity because all land–atmosphere exchanges would be computed for each patch individually but would eventually become prohibitive in computing time due to the ever-increasing number of tiles. The second approach consisted of immediately merging any newly created patch with the rest of the grid cell, thereby leading to a single tile. This Average approach would have the smallest computing burden but would entirely neglect subgrid cell dynamic heterogeneity.
The third approach consisted of keeping track indefinitely of all patches created by stand-clearing disturbance events, but grouping these patches into a much smaller number of tiles for the evaluation of land–atmosphere exchanges. In this Landscape Fractions approach, the patches having a similar (i.e., falling between fixed and predefined thresholds) were temporarily grouped into a given landscape fraction tile (henceforth LF, in roman). Land–atmosphere exchanges for each LF were then computed based on the mean area-weighted state of all patches composing the tile. Figure 1 illustrates the procedure for a grid cell with 10 patches: nine patches created by fire events (the ellipses) and one patch for the yet undisturbed part of the grid cell. Here is computed independently for each patch, along with the resulting value for . The 10 patches are then sorted by order of increasing . The patches are finally assigned to the different LF, with each LF covering a given range of values (LF1 necessarily starts at , whereas the last LF necessarily ends at ). The sorting and assignment process is repeated periodically, so that a typical patch will progressively go from LF1 to LF2 (and so forth) over time. The maximum number of landscape fractions in a given grid cell (henceforth , in italics) is also fixed and predefined.
The Landscape Fractions would be intermediate between the Full and Average approaches, both in precision (as illustrated in section 3 with the HLM) and in computing time, because land–atmosphere exchanges would be computed independently for tiles instead of many more tiles (one per patch) for the Full approach or a single homogeneous tile for the Average approach. Computing these exchanges independently for different PFT is now common in climate models, but each PFT is assumed homogeneous within its own tile (Bonan et al. 2002; Best et al. 2011; Lawrence et al. 2012; Li and Arora 2012; Reick et al. 2013). The computation of at a finer level (i.e., individual patch) than the computation of all land–atmosphere exchanges (i.e., an LF tile made up of many patches) is similar to the scaling down philosophy of the ORCHIDEE-FM (Bellassen et al. 2010) and CABLE-POP (Haverd et al. 2013, 2014) models, except that these previous models consider a single tile per grid cell for the computation of land–atmosphere exchanges.
In the context of albedo computation with the HLM, the three approaches differed only in their evaluation of [Equation (1)], which reflected the application of the versus function [Equation (2)] to each patch (Full), to a single grid-mean amount of stem carbon (Average) or to each LF (Landscape Fractions). Mathematical details on the computation of and the dynamics are provided in supplemental material 1.1. Since we assumed that the forest was in equilibrium with the fire or logging disturbance regime and that the total stem carbon in the grid cell was the same for the three approaches, the HLM results for the Landscape Fractions approach would be equal to 1) the Full approach if were equal to the number of patches or 2) the Average approach if . We address in section 3 the choice of the specific thresholds between consecutive LF (see Figure 1). The only difficulty in applying the HLM resided in the formally infinite number of patches for the equilibrium fire regime; we explain in supplemental material 1.2 how we circumvented this issue.
2.1.4. Input data
Computing RSOL in the HLM [Equation (1)] required different sets of input data on , , and . Besides the approach chosen, the value of depended upon (which determined the age-class distribution of patches) and . This last function involved , which was itself computed from the values of , a, and f. Except for the function, we employed empirical measurements from Canadian boreal forests or values from calibrated terrestrial carbon models to generate the HLM input data (see supplemental material 1.3 for more details). For a, we used a constant value of 0.30, which seemed appropriate for black spruce in boreal Canada (Gower et al. 2001). For other input data, we identified representative values as well as upper and lower values for sensitivity assessments.
For the turnover-related parameter f, we used values of 0.980, 0.985, and 0.990. required an input average value over the stand life cycle as well as different parameters for the shape of the as a function of age. We used values of 0.15, 0.25, and 0.35 kgC m−2 yr−1 for . The different shapes for are illustrated in Figure 2, along with the daily input values for , , and . Note that the functional form we used for allowed us to solve analytically the dynamics for each patch (see supplemental material 1.1.2).
The lower, representative, and upper values we used for were 0.5%, %, and 1.0%, respectively. Finally, for the upper and lower values of , we used linear functions enveloping the representative function given by Equation (2) and kept the same maximum limit of 0.975. The function used for the upper values had the same slope as Equation (2) for , whereas the function used for the lower values went from zero to the point where the two terms on the right-hand side of Equation (2) are equal (see Figure 3).
2.2. List of simulations
For each disturbance type (logging or fire), the total number of combinations of the different sets of input data is (three different sets of values for each of the seven variables). We therefore decided to divide the input variables into two groups and apply a two-level simulation design to assess the sensitivity of the HLM results (see Table 1).
The three level-1 variables (, the covarying albedo values, and ) would likely show the most substantial year-to-year variability in real climate models. We assessed the sensitivity of HLM results for all their possible combinations, thus leading to different simulations. For all 27 combinations, the level-2 variables kept their representative values. The combination of level-1 input variables that resulted in the largest (smallest) difference in RSOL between the Full and Average approaches is henceforth designated as the Maximum (Minimum) case, whereas the combination for which all level-1 input variables took their representative values is designated as the Standard case. For the four level-2 variables [sol, f, and the shape of the and functions], we assessed the sensitivity of HLM results by varying only one input variable at a time, the other three keeping their representative values. This led to 2 (upper or lower) 4 (input variables) 8 combinations. For each combination, we selected the Standard, Maximum, and Minimum cases for level-1 input variables, thus leading to different simulations.
For each disturbance type, we performed the 51 simulations (27 for level-1 variables and 24 for level-2 variables) for 1) the Full approach, 2) the Average approach, and 3) the Landscape Fractions approach using equal to 2, 3, 4, and 5.
3.1. HLM results
For the same level of disturbance , logging led to a younger set of patches than fire. Given that the albedo of a patch decreased as increased and that generally (but not always, due to decreasing NPP for older stands) increased with age, RSOL was usually higher for logging than for fire, regardless of the approach (Full, Average, or Landscape Fractions). However, these differences in RSOL between fire and logging were relatively small (i.e., always W m−2, with mean values W m−2) when comparing results for the same combination of input data and the same subgrid approach. Consequently, we present in the main text the HLM results for logging only and refer readers to supplemental material 2.1 for a comparison with fire results.
Henceforth, stands for RSOL evaluated through approach X minus RSOL evaluated through approach Y. The across the 51 combinations of input data varied from 0.0 to 8.0 W m−2, with an average of ~3 W m−2 (mean of 3.3 and median of 3.1 W m−2). Both extreme values came from the use of the lower function (Figure 3), which was used for the purpose of sensitivity assessment and does not appear realistic (see supplemental material 2.2 for additional analyses of the HLM results). For the majority of cases (47 out of 51), was between 1.5 and 6.2 W m−2. The fact that all values were positive means that the Average approach always underestimated the reflection of solar radiation. This systematic underestimation came from the downward concavity of the function, which in turn came from the fact that the tree cover over a given area must level off for high amounts of tree carbon. This bias of ~3 W m−2 corresponded to ~15% of RSOL computed according to the Full approach. The magnitude of thus means that resorting to the Average approach to account for subgrid dynamic heterogeneity can lead to substantial inaccuracies in energy budget simulation, at least over boreal forests.
Figure 4 shows as a function of the number of used in the Landscape Fractions approach. For , the Landscape Fractions approach reduced to the Average approach. For , the results depended upon the range of values included in each LF. For each of the 51 simulations performed with each value of (greater than 1), we tried all the possible thresholds to separate consecutive LF and selected, for each instance, the values for which was minimized. In other words, Figure 4 shows the results that corresponded to the best outcomes from the Landscape Fractions approach with regards to the choice of the thresholds. The results for different values of came from the use of a piecewise linear function. For the more realistic combinations of input data, decreased quickly as increased. For , was always less than 0.6 W m−2 and had a mean value smaller than 0.3 W m−2. Using , the error was always less than 0.2 W m−2 and had a mean value smaller than 0.1 W m−2. Figure 4 also shows the specific results for the most representative combination of input data (the symbol), that is, the combination for which each set of input data took its representative values.
The previous results illustrate the potential of the Landscape Fractions approach to substantially improve the RSOL results compared with the Average approach, even for a small value. Given that the previous results corresponded to the optimal outcomes from the Landscape Fractions approach, one may, however, wonder whether the improvement persists when using nonoptimal values for thresholds. Figure 5 suggests that this is indeed the case. The figure shows , with , as a function of the values of the thresholds (the fact that the second threshold is necessarily greater than the first threshold explains the triangular shape of the results domain). Note that the results were for the most representative combination of input data only. For the optimal thresholds (shown with a symbol in Figure 5), was equal to 0.31 W m−2. For the worst possible choices of thresholds, was ~3 W m−2 greater than for the optimal thresholds. Over more than half of the thresholds domain, however, W m−2 (see the zone delimited by the dashed line).
The robustness of the improvement brought by the Landscape Fractions approach was valid not only for the most representative combination of input data but for all 51 simulations. Furthermore, the domains (in terms of threshold values) of substantial improvement overlapped across all simulations. As long as reasonable values were chosen, the exact same set of nonoptimal thresholds therefore led to improvements that were robust to any change in input data. We tested two such sets and found that for was always smaller than 0.8 W m−2, with a mean value of 0.3 W m−2 across the 51 combinations of input data (see supplemental material 2.3). Finally, we also validated that the robustness (and not only the magnitude) of the improvement increased with .
3.2. Limitations of the HLM
Our computation of through Equation (1) implicitly assumed that and did not vary with the age of each patch. This assumption was a reasonable first-order approximation for treeless areas, especially for the snow-cover period, but not for tree-covered areas (Amiro et al. 2006). Accounting for the age dependence of albedo in tree-covered areas should increase (see supplemental material 2.4.1) and have minor impacts on the other main HLM outcomes.
Our model was dynamic at the scale of individual patches, but we performed the computations for a grid cell in equilibrium with constant disturbance regimes. Fire, however, does fluctuate markedly over long (Bergeron et al. 2010) and short (Stocks et al. 2002) time scales. The strict enforcement of even-aged forest management over various decades also seems highly unlikely. Nonconstant disturbance regimes should thus prevent any grid cell from ever reaching a perfect equilibrium. Nevertheless, the out-of-equilibrium grid cell can be conceived of as being in a transition between conditions similar to the ones we modeled. We therefore consider that the main HLM outcomes likely apply to out-of-equilibrium conditions, although the specific numerical results would certainly change.
Except for the spatial heterogeneity caused by logging or fire, we assumed that all factors were uniform in space as well as constant or periodic over time. The main reasons for this simplification were to keep results intelligible and be able to solve analytically the carbon dynamics of the HLM. Yet even in the state-of-the-art climate models for which our study is relevant, many abiotic (e.g., soil properties) and biotic (e.g., traits variation among and within tree species) terrestrial heterogeneities remain underrepresented (Moorcroft 2006).
Instead of simply using the most representative values, we performed a range of simulations to show that the main HLM outcomes were robust to changes in the specific values of input variables. Nonetheless, the choice of these values still involved some uncertainty. In the supplemental material, we further discuss the uncertainty related to the function (supplemental material 2.4.2) as well as the additional analyses we performed (supplemental material 2.4.3) to show that the main HLM outcomes were robust to reasonable variations in the other input data.
The HLM results presented previously led to the following outcomes. 1) On average (both mean and median values), the difference between fully accounting for the subgrid dynamic heterogeneity and performing computations based on a single average value was ~3 W m−2, reaching values higher than 5 W m−2 in 8 out of the 51 comparisons we performed. 2) The intermediate Landscape Fractions approach considerably improved the results compared with the Average approach, with an accuracy that increased quickly with the user-selected value. 3) This improvement was robust to the use of nonoptimal (but reasonable) threshold values to determine the domain of each LF. 4) When the disturbance rate was the same, between approaches were quantitatively similar, and qualitatively equivalent, for even-aged logging management and spatially random fire events of constant size.
4.1. Consequences of subgrid cell dynamic heterogeneity
4.1.1. Net radiative forcing from stand-clearing disturbances
The HLM results showed that the approach used to account for subgrid dynamic heterogeneity can have a substantial impact on the estimation of over boreal forests. Previous studies found that large changes in tree cover over boreal forests lead to an albedo-induced cooling similar to, or greater than, the warming effect from increased atmospheric carbon (Betts 2000; Claussen et al. 2001; Randerson et al. 2006; Bala et al. 2007; Bathiany et al. 2010; Arora and Montenegro 2011; Pongratz et al. 2011; O’Halloran et al. 2012). Not accounting adequately in climate models for the albedo increase that follows the temporary loss of tree cover could thus give misleading results on the net climatic impact (i.e., warming or cooling) of stand-clearing disturbances, as illustrated by the additional analysis we present here.
To estimate the magnitude of this potential error, we compared the radiative forcing of carbon emission and albedo change resulting from fire simulated in the HLM (see supplemental material 3.1 for more details). We assumed that the entire region composed of boreal ecozones in Canada was either completely free from any disturbance (the Old Growth scenario) or in equilibrium with a spatially random fire disturbance regime having equal to the 1959–97 average burn rate of the region (the Disturbed scenario). We then used the HLM to quantify the changes in atmospheric carbon content and between the two scenarios. We found that the net radiative forcing caused by fire (i.e., from the changes between the Disturbed and Old Growth scenarios) was negative for both the Full and Landscape Fractions approaches but positive for the Average approach (see Table 2). This qualitative difference came only from the major underestimation of by the Average approach because the HLM assumed equal carbon emissions for the three approaches.
Simulating boreal forests stand-clearing disturbances in climate models through the Average approach could thus bring a bias toward a warming impact. Although the estimation of radiative forcing with the HLM was relatively crude, we consider that the results obtained would likely have the same order of magnitude with state-of-the-art methods. Indeed, the results were similar to the ones of Randerson et al. (2006) when applying the HLM to the specific fire event they studied (see supplemental material 3.2).
4.1.2. Additional impacts
The albedo estimated for snow-covered conditions varies considerably among climate models, which strongly affects the simulated strength of the snow-albedo feedback for current conditions and future projections (Qu and Hall 2014). Despite having consequential impacts on future warming in the Northern Hemisphere, this source of uncertainty has apparently not been reduced over the last decade (Qu and Hall 2014). Essery (2013) concluded that the main cause of this uncertainty does not reside in the mathematical form of the vegetation–snow albedo algorithms but rather in the simulated distribution of vegetation cover and/or the specific values of albedo parameters. The inadequate simulation of the south-to-north decline in tree cover is an example of such shortcomings in current climate models (Loranty et al. 2014). Our results suggest that even if the vegetation distribution and the albedo parameters were appropriate, climate models would need to properly account for the subgrid heterogeneity created by stand-clearing disturbances in order to precisely estimate albedo, at least over disturbance-prone portions of boreal forests.
We suggest that the carbon cycle could also be affected by the approach taken to implement stand-clearing disturbances in climate models. Using the CTEM DGVM coupled to the CLASS land surface scheme (Arora and Boer 2005), Li and Arora (2012) examined the differences in simulated energy and carbon exchanges between the mosaic and composite approaches under the static coexistence of two PFT; for the simulations performed in boreal or temperate locations, the grid-level differences in the major energy fluxes between the two approaches were all below 10%, yet the differences in the equilibrium values of carbon-related variables (, biomass, and soil carbon) were mostly ~20%–40%. These major changes arose from the minor differences in the physical environment (mainly net radiation, temperature, and moisture) experienced by the PFT between the two modeling approaches. In the context of stand-clearing disturbances, the previous results suggest that limited differences in biogeophysical fluxes like can translate into consequential impacts on the simulated carbon dynamics during forest regrowth.
4.2. Implications for current climate models
4.2.1. Representing subgrid cell dynamic heterogeneity
Many models referred to as DGVM actually perform only a subset of all the computations required for full integration within climate models (Prentice et al. 2007; Quillet et al. 2010). For the sake of clarity, we henceforth designate as Dynamic Vegetation–Land Surface Models (DVLSM) the process-based models that compute the land–atmosphere exchanges of energy, carbon, water, and momentum required by climate models, while allowing vegetation to change dynamically (both the spatial distribution and the seasonal cycle) as a function of the climatic conditions. In other words, DVLSM correspond to the union of first- to third-generation land surface schemes, as defined by Pitman (2003), with dynamic vegetation. To the best of our knowledge, only one Earth system climate model currently uses a DVLSM (i.e., SEIB–DGVM) that includes a second-generation treatment of vegetation interactions by allowing many individuals of the same PFT in a single grid cell. We therefore focus on DVLSM using a first-generation representation of vegetation before coming back to this exception below.
Few studies in Earth system climate models have directly addressed the impacts of stand-clearing disturbances. Lawrence et al. (2012) evaluated the climatic impacts of logging (along with long-term changes in land cover) at the global scale, from 1850 to 2100. They represented subgrid dynamic heterogeneity through the Average approach within a PFT-tiling nested framework, resorting to MODIS and potential vegetation data to derive transient PFT coverage datasets. Landry et al. (2015) estimated the effect of different fire scenarios on global carbon stocks and the global atmospheric surface temperature over the 2015–2300 period. They once again represented subgrid dynamic heterogeneity through the Average approach, within a PFT-tiling framework based on simulated competition. According to the HLM results, these two studies may thus have underestimated the albedo-induced cooling caused by stand-clearing disturbances, at least over boreal forests. As far as we can tell, all first-generation DVLSM representing fire and/or logging in Earth system climate models (e.g., ORCHIDEE, JSBACH, and CLM4) resort to the Average approach, except LM3V. Unfortunately, the Average approach misrepresents the transient biogeophysical impacts resulting from temporary loss of forest cover as we showed for with the HLM because the disturbed patch is averaged out with mature forest. This inadequate representation of subgrid dynamic heterogeneity also likely affects other biogeophysical exchanges and, as explained above, possibly the carbon cycle itself. As for LM3V, it creates a new “secondary” tile for logging each year but treats fire through the Average approach (Shevliakova et al. 2009). Moreover, LM3V follows a maximum of 12 secondary tiles per grid cell, which is probably insufficient to adequately represent subgrid heterogeneity in regions where logging happens each year and vegetation grows slowly.
As mentioned above, SEIB–DGVM is, to our knowledge, the only second-generation DVLSM that has been included in an Earth system climate model (Watanabe et al. 2011). As many stochastic gap models, SEIB–DGVM simulates the coexistence of various PFT in a single tile of ~0.1 ha (Sato et al. 2007). Accurately representing within such a small area undisturbed and various burned or logged patches does not seem possible (e.g., a few surviving trees could shield most of the tile’s ground). Due to the stochastic nature of modeled processes, gap models need to simulate various replicates, typically 10 or more, for each grid cell. Some of these replicates could be entirely burned or logged in order to accurately represent the impacts from stand-clearing disturbances. This would however require increasing the number of replicates following each disturbance event, thereby increasing the computing time that is already an issue with gap models (Scherstjanoi et al. 2013). Consequently, accounting for the whole age-class distribution of patches created by stand-clearing disturbances seems challenging in SEIB–DGVM and similar DVLSM, especially in slow-growing boreal forests where decades can elapse before canopy closure. Incidentally, we note that the SEIB–DGVM built-in fire module (Sato et al. 2007) was not activated for the CMIP5 simulations (Arora et al. 2013). We also note that the ED model, which now computes for each patch the land–atmosphere fluxes required from DVLSM (Medvigy et al. 2009), basically implements the Full approach by default, except that patches are merged together to minimize computing time when they become sufficiently similar (Fisher et al. 2010). However, the ED model has apparently not been integrated yet into an Earth system climate model or even used for global simulations as a stand-alone DVLSM, likely reflecting its computing demand.
Finally, increasing the spatial resolution of DVLSM would generally be insufficient to resolve subgrid dynamic heterogeneity in global simulations; for example, individual fire events of less than 500 km2, each covering less than 7% of a 1° × 1° grid cell (centered on 50°N), are responsible for 55% of the mean area burned in the boreal ecozones of Canada (Stocks et al. 2002).
4.2.2. Implementing the Landscape Fractions approach
Given the potential biases of the Average approach and the high computing time of the Full approach, our HLM results suggest that approaches similar to the Landscape Fractions could constitute an interesting middle point for the representation of stand-clearing disturbances in first-generation DVLSM. Assigning patches to specific LF and updating their status each month or year should involve a small computing burden so that the additional time required by the Landscape Fractions approach should increase roughly linearly with . The exact rate of increase would depend upon how much time the specific DVLSM spends computing the frequent land–atmosphere exchanges. In the case of IBIS, this fraction can be around 60% for a standard simulation. The accuracy of the results, on the other hand, seems to increase exponentially with (see Figure 4). Based on the HLM, would appear as a good trade-off between accuracy and computing time. For this choice, the total computing time in IBIS would approximately double [3 () 60% (fraction of the computing time spent on land–atmosphere exchanges) 40% (rest of the computing time) 220%]. Please refer to supplemental material 3.3 for more information on the implementation of the Landscape Fractions approach in first-generation DVLSM.
Stand-clearing disturbances lead to the temporary removal of tree cover over extensive areas each year but are then followed by forest regrowth. Given that these events are only starting to be the subject of dedicated studies in Earth system climate models, the timing seems appropriate for a thorough consideration of the consequences from different approaches that could be implemented to account for the resulting subgrid cell dynamic heterogeneity.
In this study, we specifically assessed three such approaches. The Full approach, which consisted of creating an individual tile for each new patch created by a disturbance event, would result in the most accurate representation of subgrid dynamic heterogeneity but would quickly become prohibitive for long simulations. The Average approach, which consisted of immediately merging any new patch with the rest of the (forested part of the) grid cell, would be a simple and time-efficient strategy but would actually neglect subgrid dynamic heterogeneity. The third approach that we assessed aimed to be intermediate both in accuracy and computing time. The basic idea behind this Landscape Fractions approach was to temporarily group the patches having a similar proportion of tree cover into the same landscape fraction (LF) tile. In climate models, the status of individual patches would then be updated infrequently (monthly to yearly time scales), whereas the time-consuming computations of land–atmosphere exchanges of energy, carbon, water, and momentum performed at a hourly time scale would be done in a much smaller number of LF tiles.
We illustrated the consequences of the three approaches in the Heterogeneous Landscape Model (HLM) we developed for this study by computing the annual reflection of solar radiation () over a boreal evergreen forest in equilibrium with a constant fire or logging disturbance regime. Even if we assumed the same total amount of carbon for the three approaches, we found that the difference between the Full and Average approaches was ~3 W m−2 on average (i.e., about 15% of ), often reaching values higher than 5 W m−2. The results from the Landscape Fractions approach converged quickly toward the results from the Full approach as the number of LF increased and were robust to the use of nonoptimal values for the thresholds separating consecutive LF. Simulations with the HLM therefore brought support to the idea that the Landscape Fractions approach could lead to substantial improvements in accuracy over the Average approach, at a fraction of the computing time required by the Full approach, thereby constituting an intermediate strategy between the simplest and most accurate approaches, respectively. We also found that, overall, the difference in the age-class distribution of patches between the fire and logging disturbance regimes had a negligible impact on the values between the approaches. Finally, we performed an additional analysis with the HLM to assess the net radiative forcing resulting from fire across the boreal ecozones of Canada, as a function of the approach taken to represent subgrid dynamic heterogeneity. The outcome suggested that simulating boreal forest stand-clearing disturbances in coupled climate–terrestrial vegetation models through the Average approach would bring a considerable bias toward a warming impact.
We do not believe that precisely accounting for subgrid cell dynamic heterogeneity should be a priority under all circumstances. When the objective is to compare widely diverging scenarios of anthropogenic emissions, the biases resulting from the Average approach are likely acceptable because they probably cancel each other to a substantial extent among scenarios and are anyway much smaller than other sources of uncertainty, for example, the simulation of clouds and their interactions with aerosols. However, when the objective is to directly assess stand-clearing disturbances, neglecting the resulting subgrid dynamic heterogeneity might be misleading, particularly when studying the net climatic impact from fire or logging in boreal forests.
We thank David Price for helpful discussions about the development of the HLM and for providing us with the monthly mean incoming solar radiation data. We also thank Günther Grill, Gabriela Jamett, and Bano Mehdi for their comments on a previous version of the manuscript as well as Damon Matthews and Julia Pongratz for discussions about the Landscape Fractions approach. Thoughtful comments from the two reviewers have further helped us improve the manuscript. J.-S. L. was supported by an NSERC postgraduate scholarship (CGSD) and an FRQNT doctoral scholarship (B2).
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/EI-D-15-0022.s1.
Current affiliation: Department of Geography, Planning and Environment, Concordia University, Montréal, Canada.
Current affiliation: Liu Institute for Global Issues and Institute for Resources, Environment, and Sustainability, University of British Columbia, Vancouver, Canada.
Current affiliation: Earth and Environmental Sciences and Biology, Irving K. Barber School of Arts and Sciences, University of British Columbia, Kelowna, Canada.
This article is included in the Biogeophysical Climate Impacts of Land Use and Land Cover Change (LULCC) special collection.