**1. Introduction**

Boreal forests are the northern-most forested biomes, and are expected to be sensitive to climate change and forest fire disturbances [1]. By 2100, the climate in Northern high latitudes (including North America and Eurasia) is expected to increase by 1.4 to 5.8 ◦C, about 5 times greater than the global mean temperature increase [2,3]. Recent studies have projected that climate change will affect the species distribution and ecosystem functions (e.g., carbon fixation) within boreal forests [4]. For instance, climate change can alter interspecific competition and tree species migration [5], which could further affect forest composition and distribution [6]. Species response to climate change varied among and within populations of different trees in the boreal forests [7]. Changes in climatic conditions (e.g., precipitation and temperature) can have direct influences on the metabolic processes, growth rates, establishment abilities, and competitive ability of trees, thus affecting overall biomass accumulation in forests [5,8]. These resulting changes in forest structure and function are expected to affect the recovery ability (ecological resilience) of boreal forests at landscape scales [9]. Additionally, climate change indirectly impacts forest traits (e.g., forest structure, composition, and ecological resilience) through its effects on fires regimes [10]. In boreal forests, climate-induced fires are frequent and widespread, and become a major factor that can distinctly affect forest successional dynamics, composition, and structure [11,12]. Johnstone et al. [11] showed that forest fires with increased severity may promote shifts from coniferous forest to deciduous-dominated forests, and substantially change landscape dynamics and ecosystem services in boreal forests. More previous studies indicated that fires have been projected to occur more frequently, burn greater areas, and have higher intensities under altered climatic conditions [13,14]. As a result of climate change altered fire regimes, forest composition and biomass dynamics, and thus ecological resilience is expected to shift [15]. Despite the growing evidence that climate change and shifting fire regimes will alter the composition, structure and biomass of boreal forests, quantification of how these two factors will impact forest ecological resilience is still poorly known.

Ecological resilience is characterized as the capacity of a forest ecosystem to recover from disturbance and maintain a stable state, supporting the recovery of structure, composition, and function equivalent to pre-disturbance states [16]. Boreal forests were remarkably resilient to disturbances, and forest species were adapted to the current disturbance regimes with long term effects [17]. For the existence of forest ecological resilience, boreal forest ecosystems thus have the capacity to absorb a spectrum of perturbations (e.g., climate change and forest fires) and to sustain its structure and function, and to maintain the forest ecosystem in a relatively stability domain [16,17]. Climate change in the past century has caused more frequent extreme climate events, such as higher temperatures, severe, and extensive droughts [15], and also has altered forest fires regimes to varying degrees [12,18]. These changing factors (e.g., climate change and climate-induced fires) will exacerbate the loss of ecological resilience in boreal forest ecosystems under long-term exposure [15,19], and may cause a catastrophic shift in forest ecosystems that is difficult to reverse, thus posing a very serious threat to regional ecological security and forest service [20]. Therefore, understanding and quantifying ecological resilience is increasingly important for forest ecosystem management, and provides a quantitative basis for exploring the issue of maintaining and improving ecological resilience.

Harvesting and planting are major anthropogenic disturbances to boreal forests. Boreal harvesting and planting alters forest composition and structure, aboveground biomass accumulation, and ecological resilience from stand to landscape scales [21], and these effects could be aggregated under future changed climate conditions [22,23]. He et al. [22] evaluated species response to harvesting and climate-induced fire in Northern Wisconsin boreal forests, and showed that increased fire frequency can significantly alter the distribution of shade tolerant species, and indicated that harvesting accelerated the decline of Northern hardwood and boreal tree species. Gustafson et al. [23] estimated the climate effects on forest composition in the South-Central Siberian region, and indicated that the direct effects of climate change were not as important as the timber harvesting effects on local virgin forests. However, there are fewer studies exploring the effects of forest managemen<sup>t</sup> schemes (harvest and planting strategies) on the ecological resilience of boreal forests.

Climate change, fire disturbance, harvesting, and planting occur at large spatio-temporal scales, which makes evaluating their effects on ecological resilience using traditional observation experimental studies challenging [24,25]. Forest landscape models (FLMs) provide a proper scientific approach for studying these issues [26]. With FLMs, we can conduct large-scale studies in which critical model parameters could be changed to explore the complex interactive effects of these extra factors on ecological resilience [27].

The objective of this research was to investigate effects of climate change, climate-induced fire regimes, and future possible forest managemen<sup>t</sup> schemes on the ecological resilience of boreal forests in Northeastern China. Specifically, we quantified (1) individual effects and (2) interactive effects of climate change and climate-induced fire disturbance on boreal forest ecological resilience, and (3) evaluated whether future possible forest managemen<sup>t</sup> schemes could mitigate the effects of climate change and climate-induced fires on ecological resilience.

#### **2. Materials and Methods**

## *2.1. Study Area*

The study area is located in the Great Xing'an Mountains, which covers nearly 2.7 million ha (Figure 1, 51◦35 to 53◦25 N and 122◦25 to 125◦35 E). The climate conditions are characterized by terrestrial monsoons with long winters and short summers, and the mean monthly temperatures range from −28 to 20 ◦C in January and July, respectively. Precipitation mainly falls in the summer, and the mean annual value is 428 mm. The elevation ranges from 173 m to 1511 m across the landscape, and the region is covered by brown coniferous forest soils. The dominant vegetation in this area is larch (*Larix gmelinii*) forests. White birch (*Betula platyphylla*) and aspen (*Populus davidiana*) are the major broad-leaved species in this region. In addition to larch and white birch, Mongolian Scots pine (*Pinus sylvestris* var. *mongolica*), and Korean spruce (*Picea koraiensis*) are also widely distributed. Dwarf pine (*Pinus pumila*) has small species communities which can be found in high latitude regions.

Forest fire is a major disturbance in the Great Xing'an Mountains. Based on the Chinese Federal Forest Service data (website: http://www.cfsdc.org), forest fires burned 519,144 ha of the landscape during the period of 1965 to 2005 in this region. For half a century, extreme fire suppression policies have changed fire regimes in this area profoundly. Previous studies indicated that fire regimes have changed from frequent and lower intensity fires to more infrequent and high intensity fires [28]. Extensive harvesting events have affected the forest structure, composition, and natural regeneration significantly in this region. According to the forest inventory data and our field investigation, coniferous dominated forests have shifted from late-seral to mid-seral stages over the landscape. Planting occurs rarely in our study area and overall has minimal effect compared to fire and harvesting. Under the long-term effects of these two typical disturbances (fire and harvesting), the boreal forests of our study region have become more fragmented, simplified, and less resilient [29].

**Figure 1.** The location of our study area. The insert map shows the location of our study region in Northeastern China.

#### *2.2. The Indicators of Ecological Resilience*

The goal of our study was to evaluate the effects of three variables (climate change, climate-induced fire regimes and future possible forest managemen<sup>t</sup> schemes) on ecological resilience of the study region. Firstly, by using LANDIS PRO model, we simulated the effects of climate variations, fire regimes, and forest managemen<sup>t</sup> schemes on these four indicators: Aoveground biomass, growing space occupied, age cohort structure, and proportion of mid and late-seral species, which can be obtained directly by LANDIS PRO model outputs. Secondly, we estimated ecological resilience through these four indicators (weighted by a function of the variance). Thus, we can investigate the effects of climate variation, fire regimes, and forest managemen<sup>t</sup> schemes on boreal ecological resilience in our study area.

Previous studies indicated that the ecological resilience can be quantified through boreal structure, composition, and functioning [15,30]. Seidl et al. [9] and Van Mantgem et al. [31] suggested that the total C storage, the rumple index, the presence of late-seral species and the proportion of older age cohorts can be served as the indicators of ecological resilience. Based on the two previous studies and the current status of our study area, we selected aboveground biomass, growing space occupied, age cohort structure and proportion of mid- and late-seral species as the indicators of ecological resilience (Table 1). The specific contents of these selected indicators are as follows: With regard to boreal forest functioning, we mainly focused on aboveground biomass (AGB). AGB plays an important role in carbon fixation, and is a vital surrogate of forest ecosystem functioning [9]. As a surrogate of forest structure, we used the growing space occupied (GSO) as an indicator of ecological resilience. The GSO is the growing space occupied by species of a specific site, and is commonly used as the crown closure measurement. Growing space primarily reflects forest structure, and can distinguish forest successional stages over large landscapes [32]. Seidl et al. [9] selected the rumple index (RI) of canopy complexity as a surrogate of vegetation structure. The rumple index is the ratio of the canopy surface area to the projected surface ground area [33]. The RI was proposed as a powerful composite index to describe vegetation structure and distinguish different stages of forest development over large areas [34], which was similar to the contents of GSO. Thus, we used GSO as the indicator of ecological resilience in this study. The measure of GSO incorporated in LANDIS PRO is calculated by the following equations:

$$\text{GSO}\_{\text{i}} = \sum\_{\text{j=1}}^{\text{longactive}\_{\text{j}}/\text{times}} \left( \left( \frac{\text{DBH}\_{\text{i}}}{\text{10 inch}} \right)^{1.605} \times \text{NT}\_{\text{j}} \times \frac{1 \,\text{lactare}}{\text{MaxSDI}\_{\text{i}} \times \text{site\\_area}} \right)$$

$$\text{GSO} = \sum\_{\text{i=1}}^{\text{number of species}} \text{GSO}\_{\text{i}}$$

where DBHj and NTj are the mean diameter and number of trees of j th diameter class for species i in inches and stems, respectively; MaxSDI is the maximum stand density index (derived from species vital attributes). Stand density index (SDI) is a basic concept in forestry. It was first developed by Reineke in 1933 and has been widely used to characterize stand density (tree crowding). Growing Space Occupied (GSO) represents an extension of SDI to meet the needs of landscape modeling, and is also a key metric within LANDIS PRO [32].

As a surrogate of vegetation component, we selected the percentage of older age cohort individual species (ACS). Forests which contain tree species with diverse age and size structures have been observed to be relatively more resilient to climate change and fire disturbance than forests with younger, less diverse age structures [31]. ACS is defined by the ratio of old age cohorts to the total trees cohorts (trees age > 60 year for broadleaf, and >100 year for conifers). For the surrogate of composition, we selected the proportion of mid, and late-seral species (i.e., larch, Mongolian Scots pine, and Korean spruce) >5 cm in DBH (diameter at breast height) (LSS).


**Table 1.** Indicators of forest ecological resilience.

In order to facilitate the comparison and calculation among different indicators, we first used the min-max normalization method to normalize all four indicators (AGB, GSO, ACS, and LSS) by each time step at both land type and landscape scales. The data normalization process is calculated using the following Formula (1):

$$X\_i = \frac{|X\_i - X\_{\min}|}{X\_{\max} - X\_{\min}} \tag{1}$$

where *Xi* is the normalized data, which ranging from 0 to 1. *Xi* is the value of the indicator at year *i*. *Xmin* and *Xmax* represent the minimum and maximum values, respectively.

We then calculated the ecological resilience at all simulated time steps by using all four indicators (AGB, GSO, ACS, and LSS). The calculated ecological resilience is a specific number ranging from 0 to 1 that quantifies the capacity of different forest stands to recover from extra disturbance and maintain a stable status. Forest stands with high resilience values have a higher ability of recovery (more resilient than other stands). The formula for this is (2):

$$R\_i = \mathcal{W}\_1 \times \overline{\overline{AG}}\_i + \mathcal{W}\_2 \times \overline{\overline{GSO}}\_i + \mathcal{W}\_3 \times \overline{\overline{AC}}\_i + \mathcal{W}\_4 \times \overline{\overline{LSS}}\_i \tag{2}$$

where *Ri* is the ecological resilience value at year *i*, higher *Ri* means higher ecosystem recovery ability (ecological resilience). *AGBi*, *GSOi*, *ACSi*, and *LSSi* represent the normalized AGB, GSO, ACS, and LSS

values at year *i*, respectively. *Wj* (*j* = 1, 2, 3, 4) are the weight coefficients, which are calculated by using the coefficient of variation method in the following Formula (3):

$$\mathcal{W}\_{\dot{j}} = \frac{\sigma\_{\dot{j}}}{\overline{\overline{X}}\_{j} \times \sum\_{j=1}^{n} \frac{\sigma\_{\dot{j}}}{\overline{X}\_{j}}} \tag{3}$$

=

where *Wj* is the weight coefficient of indicator *j*; *σj* is the standard deviation of indicator *j*; *Xj* is the mean value of the indicator *j*.

#### *2.3. Simulation Experiments Design and Data Analysis*

We designed a factorial experiment to assess the effects of climate change, climate-induced fire regimes, and forest managemen<sup>t</sup> schemes on boreal forest ecological resilience. In this factorial experiment, we set three independent variables: Climate change (current climate and climate change), different fire regimes (current fire and climate-induced fire), and forest managemen<sup>t</sup> schemes (no treatment and different harvesting and planting strategies).

The current meteorological data were derived from the meteorological center, and monthly temperature and precipitation data were included from year 1961 to 2000. We used the data derived from five related weather stations to build regression models among spatial positions, elevations, and temperature as well as precipitation in the studied area. We then calculated the mean annual temperature and precipitation of different land types by using this regression model. We used two different levels of carbon emissions scenarios (CGCM3 B1 and UKMO-HadCM3 A2) to represent future climate change in our study. The B1 scenario represents low CO2 emissions, while A2 represents high CO2 emissions [3,35]. Based on the projected data of Hadley GCM, the mean annual temperatures and precipitations would increase linearly in year 2000-2100, and after that it would enter into a stable state [36]. The historical fire regimes for our simulations were characterized by the Chinese Federal Forest Service database from 1965 to 2005. According to previous study, fire occurrences in our study region under the B1 and A2 scenarios (projected by the Hadley GCM) would increase by 30% and 200% compared to historical fire regimes, respectively [3].

We used recent harvest trends in our study area to construct the current harvest regime. To examine the effects of the current harvest regime and future possible forest managemen<sup>t</sup> schemes on forest ecological resilience, we designed eight harvesting and planting scenarios (Table 2). These scenarios include a combination of designated harvest intensity and increasing percentages of individual trees planted to the current intensity (P0) to 10% (P10), 20% (P20), 30% (P30), 40% (P40), and 50% (P50) of the mean stand density.


**Table 2.** The scenarios for different harvesting and planting strategies (HP) simulated by LANDIS PRO.

Specifically, we designed five simulated scenarios: (1) Current fire only scenario (CF1: the reference scenario, fire and succession were simulated with current fire occurrence); (2) Increased fire only scenario (CF2: compared to current fire regime, fire occurrence increased by 30%; CF3: fire occurrence increased by 200%); (3) climate change only scenario (B1F1 and A2F1: climate change and current fire regimes were simulated); (4) climate-induced fire scenario (B1F2: B1 climate and fire occurrence increased by 30%; A2F3: A2 climate and fire occurrence increased by 200%); and (5) climate-fire- forest managemen<sup>t</sup> schemes (B1F2HP: B1 climate, fire occurrence increased by 30% and harvesting and planting; A2F3HP: A2 climate, fire occurrence increased by 200% and harvesting and planting). We used a FLM to simulate 5 tree species (Table S1) at 5-year time step from year 2000 to 2300 with five replicates to reduce the stochasticity.

To examine the effects of extra factors on boreal forests, we compared the response variable, ecological resilience, under the reference scenario to the scenarios climate change only, increased fire only, and climate induced-fire, and climate-fire-forest managemen<sup>t</sup> schemes for short- (0–50 year), medium- (50–150 year), and long-term (150–300 year) simulation periods using the mean comparison method. An analysis of variance (ANOVA) was used to test the differences between the reference scenario and all other scenarios. We used the Tukey's Honestly Significant Difference (HSD) method for post-hoc analyses at all simulated periods. To evaluate the increased fire effects on boreal ecological resilience, we tested the response variable among different fire regimes at all three simulated periods. All statistical analyses were conducted using SPSS 23.0 software.

#### *2.4. Simulating Ecological Resilience from Climate Change and Disturbance*

We employed a forest landscape model to simulate forest succession under different climate, fire regimes and forest managemen<sup>t</sup> schemes, and to evaluate ecological resilience of boreal forests (Figure 2). LANDIS PRO is a spatially explicit landscape model, and can be used to simulate forest dynamic over large spatial and temporal scales with user defined resolutions (10–500 m) [34]. LANDIS PRO records density and size information for each age cohort by species within raster cells enabling the model to directly output spatially explicit stand information (e.g., density, basal area, and aboveground biomass). This data structure enables forest inventory data to be directly used for model initialization and parameterization. LANDIS PRO can simulate tree growth, species establishment, mortality, and species resources competition within each raster cell, and also simulate seed dispersal, forest management, and natural disturbance across the whole landscape. LANDIS PRO stratifies the entire landscape into relatively homogenous land type units based on climate, soil, terrain, and other environmental factors. Species establishment probability (SEP) is a key input parameter of LANDIS PRO. SEPs are obtained based on responses of each species to specific microenvironment factors such as soil moisture, soil N, soil C, and local climate. LANDIS PRO uses SEPs as inputs to indirectly capture the spatial variability of climate. Species with high SEPs have a higher probability of establishment. The SEPs of specific species are derived from previous LANDIS modeling studies or a gap model (e.g., LINKAGES).

**Figure 2.** The framework for model-coupling and sub-methods used to evaluate ecological resilience. SEPs: species establishment probability, AGB: aboveground biomass, GSO: the growing space occupied, ACS: age cohort structure, and LSS: proportion of mid and late-seral species.

In the fire module, fire disturbances are simulated based on specific input parameters (e.g., mean fire return interval, mean fire sizes) [37]. Fire disturbance is simulated within spatially defined fire

regime units which have parameterized ignition rates, fire size distributions, and mean return intervals. The fire module includes three major components which are fire occurrence, fire spread, and fire effect simulation. Forest fires effects are characterized by bottom-up disturbance, and younger trees are more susceptible to fire than older ones in model simulation.

In the harvest module, forest harvesting is simulated using a managemen<sup>t</sup> area map and forest stand map. The managemen<sup>t</sup> area map and the stand map provide boundaries for specific harvest events to occur. Clear-cutting, thinning (from above or below), and group selection harvesting can be specified to execute harvest events in the harvest module [38]. By using those parameters, many common forest managemen<sup>t</sup> schemes can be simulated in this module.

#### 2.4.1. LANDIS Model Initialization and Parameterization

Parameterization LANDIS PRO included two aspects: Non-spatial parameters (species' vital attributes, SEPs, site-level fire disturbance, harvest scenarios, and planting parameters), and spatial parameters (species composition map, land type map, managemen<sup>t</sup> area map, stand map, and fire regime unit map). Five of the most common tree species were simulated in LANDIS, which account for more than 90% of the total forested land [29]. Species' vital attributes were derived from previous studies in the same or similar regions, field investigation, and consultation with local experts [39,40]. We used land use data, Landsat imagery, slope, and aspect maps to classify the study area landscape into six land types. We then used LINKAGES to simulate the response of forest species under current and climate change scenarios within each land type, and used the simulated individual species biomass to estimate SEPs for each simulated species. We modeled the SEPs of five tree species under different climate conditions by each land type (Table 3). The initial SEPs were estimated by current climate (1961–2000), and the SEPs for future scenarios were projected by climate change data (2010–2099). The SEPs were assumed to change linearly in 2000–2100, and held constant after 2100. Specifically, we used LINKAGES to simulate the individual biomass of different tree species to both current and climate warming scenarios within each land type. The individual species biomass was used to estimate the SEPs for specific species. The SEPs are calculated by the following equation:

$$SEP\_{ij} = \frac{b\_{ij} \cdot b\_{ij}'}{\max\left\{\sqrt{\sum\_{j=1}^{n} b\_{ij}}, \sqrt{\sum\_{j=1}^{n} b\_{ij}'}\right\}}$$

where *bij* and *bij* are the biomass of species *i* on land type *j* under current and warmer climate, respectively, *SEPij* is the species establishment probability of species *i* on land type *j* under current and warming climate [36,41].

The forest composition map was obtained based on forest stand maps, forest inventory data, and field data, which included species spatial location, number of trees and age cohort information. The forest stand map was acquired in the 2000's was a GIS based file, that provided stand site boundaries, species composition, structure, and the average age of the specific polygons. We derived sample plots investigated during the 2000's, which provided number of trees and age cohorts by species from the China National Forest Inventory Second and Third Tier data. We integrated the forest stand map (vector format) and forest inventory data (stand information) to derive the initial forest composition map. To reduce computational resources, all those input maps needed by LANDIS model were converted to a resolution of 90 m × 90 m cell size (2217 columns × 2609 rows) by using ESRI ArcGIS software.

Fire regime parameters and a fire regime unit map were required for the fire module. The fire regime unit map was used to identify areas with heterogeneous fire properties across the landscape, and fire characteristics in this region were mostly related to soil moisture, terrain, climate, and vegetation traits, which were closely related to the classification of land types, and thus we used the land type map as the fire regime unit map in our study area. The current fire occurrence for our simulations was parameterized based on data from the historical fire database recorded from 1965 to

2005. Based on the database, we calculated the current fire regime parameters (e.g., return interval, ignition rate, and mean fire size) of each fire regime unit. The future fire regimes were characterized by changing fire occurrences under different climate scenarios based on previous work [3]. The boreal forests in our study area have been exploited since the 1950's, and timber harvesting has extensively altered forest composition, structure, and age cohort. Consequently, to maintain forest ecosystem function and sustainability, timber harvesting has been restricted by a natural forest conservation project since 1999. Mongolian Scots pine and Korean spruce were extensively cut because of their high economic value and stands typically reestablished with larch. At present the local forestry bureaus have attempted to actively protect the remaining stock of these two species. In accordance with current harvest policy, the harvested species were larch, birch, and aspen, whereas pine and spruce were not harvested. The predominant harvest type in our study area was thinning from below, and all harvest scenarios were processed by removing the smallest trees first. We simulated the current harvest activities by using a basal area controlled harvest method (tree species were removed from a stand until a specific target basal area value was reached) followed by planting in permitted areas.


**Table 3.** SEPs for each available land type under current climate and climate change scenarios.

1 C: current climate condition; B1 and A2: climate change scenarios.

#### 2.4.2. Model Calibration and Verification

Simulated results (e.g., species composition, tree density, basal area, and aboveground biomass by species for each cell and time-step) from LANDIS PRO can be directly compared with forest inventory data as a method of calibration and validation [32]. In order to parameterize the initial forest landscape accurately, we used 70% of the inventory plots (investigated in 2000s, consists of the number of all trees and age cohorts) and the stand map (a GIS file) to initialize the forest composition map at year 2000, and then simulated the model for ten years. We iteratively adjusted species' growth curve (an essential input parameter used to control tree growth and calculate species biomass) to make the initialized forest stand information match the remaining 30% of the forest inventory data at year 2000. We then calibrated the number of potential established seeds (a parameter related to tree density and basal area) until the simulated results for year 2010 was similar to field data for the year 2010. This calibration ensured that species' growth curves and the number of potential established seeds was suitable for our study area [39]. To evaluate the simulated landscape at year 2010, we used a scatter plot of the observed density and basal area vs. the simulated density and basal area. We first selected 322 raster cells from the simulated landscape at year 2010, and then the density and basal area were extracted from selected cells to compare with forest inventory data. Likewise, the forest inventory data (322 plots, investigated in 2010s) were also converted to the total density and basal area.

To verify simulated fire on the forested landscape, we compared the model results with field data at different simulated periods. We ran the current fire only scenario (CF1) for 300 years, and randomly

selected 40 fires with low intensities (more than 90% fires occurred at this level in our study area) from different years and locations from the LANDIS PRO output. We then inventoried 40 field sites (8 field sites per each age group) that were actually burned 5, 10, 15, 20, and 25 before the year they were sampled. We set five plots with 20 m × 20 m in each field site. We then recorded all the individual trees with basal diameter above 1 cm level in each plot. Tree number and DBH by species were measured at each plot, and these plot data were converted to density (trees/ha) and basal area (m2/ha). We statistically compared tree density and basal area of the 40 simulated fires with 40 corresponding fires sampled in the field, respectively.

To ensure the simulated results more authentic, we also compared the simulated aboveground biomass to previous studies, which conducted plot surveys in similar region at landscape scales. We used the currently available data for model evaluation. While predicted results under climate change scenarios over next 280 years cannot be verified by filed inventory data, simulated successional and stand dynamic trends have been confirmed by other studies conducted at similar regions [36,42].
