Predicting Potential Fire Severity Using Vegetation , Topography and Surface Moisture Availability in a Eurasian Boreal Forest Landscape

Severity of wildfires is a critical component of the fire regime and plays an important role in determining forest ecosystem response to fire disturbance. Predicting spatial distribution of potential fire severity can be valuable in guiding fire and fuel management planning. Spatial controls on fire severity patterns have attracted growing interest, but few studies have attempted to predict potential fire severity in fire-prone Eurasian boreal forests. Furthermore, the influences of fire weather variation on spatial heterogeneity of fire severity remain poorly understood at fine scales. We assessed the relative importance and influence of pre-fire vegetation, topography, and surface moisture availability (SMA) on fire severity in 21 lightning-ignited fires occurring in two different fire years (3 fires in 2000, 18 fires in 2010) of the Great Xing’an Mountains with an ensemble modeling approach of boosted regression tree (BRT). SMA was derived from 8-day moderate resolution imaging spectroradiometer (MODIS) evapotranspiration products. We predicted the potential distribution of fire severity in two fire years and evaluated the prediction accuracies. BRT modeling revealed that vegetation, topography, and SMA explained more than 70% of variations in fire severity (mean 83.0% for 2000, mean 73.8% for 2010). Our analysis showed that evergreen coniferous forests were more likely to experience higher severity fires than the dominant deciduous larch forests of this region, and deciduous broadleaf forests and shrublands usually burned at a significantly lower fire severity. High-severity fires tended to occur in gentle and well-drained slopes at high altitudes, especially those with north-facing aspects. SMA exhibited notable and consistent negative association with severity. Predicted fire severity from our model exhibited strong agreement with the observed fire severity (mean r2 = 0.795 for 2000, 0.618 for 2010). Our results verified that spatial variation of fire severity within a burned patch is predictable at the landscape scale, and the prediction of potential fire severity could be improved by incorporating remotely sensed biophysical variables related to weather conditions.


Introduction
Wildfires, ignited by human or natural agents, are crucial disturbances in the boreal forests of Eurasia and North America [1][2][3].Wildfires can strongly influence regional land surface processes such as carbon cycling [4,5] and energy and water budgets [6,7].Severe burns can result in tree mortality and soil erosion, thereby degrading ecosystem functions [8].Nevertheless, there is a growing consensus that forest wildfires can also provide a unique opportunity for ecosystem Forests 2018, 9, 130 2 of 26 restoration [9,10].In particular, fire severity can exert profound impacts on the successional trajectories of the early post-fire vegetation [11,12], which may ultimately determine the future forest structure and function [13][14][15].Many recent studies have indicated that warming and drying climates tend to shift the current fire regimes toward more frequent, large burns with high severities [16][17][18].Understanding the causal mechanisms of fire severity patterns is essential for mitigating the adverse effects of fires, and for maintaining beneficial ecosystem functions and services.
Fire severity is generally defined as the magnitude of ecosystem changes caused by fire [19,20].It is often based on metrics obtained from the field that represent short-term fire effects in the immediate post-fire environment (e.g., tree mortality and soil organic matter loss) [19,21], or based on the validated relationships between remotely sensed spectral indices (e.g., normalized burn ratio (NBR) and differenced normalized burned ratio (dNBR)) and field-measured metrics [19][20][21][22].It has been shown that variations in fire severity patterns are at the heart of how an ecosystem will change in response to fire events [23][24][25].Spatially explicit fire severity maps can assist resource managers in evaluating post-fire biomass loss and developing sound strategies for ecological restoration [26,27].Numerous studies that incorporate field investigation and remotely sensed images have been conducted to map fire severity variability within burned patches [19,20,27,28].However, a potential fire severity map may be even more useful in allowing managers to anticipate the hotspot areas that are likely to burn severely, and thus to prioritize resources for fuel treatments for those areas [29,30].
Fire behavior is regulated by three major controlling factors (i.e., weather, topography, and fuels) that form a fire triangle at spatial scales ranging from a few hectares to thousands of square kilometers [31].Fuel composition and fuel loading interact with terrain features and fire weather to influence burning duration and heat flux, which ultimately determine the spatial pattern of fire severity [20,32].Various methods have been developed to represent the cause-and-effect relationships between environmental factors (drivers), fire behavior (process), and fire effects on an ecosystem (patterns).Those methods generally fall into two categories: physical models (e.g., FIREHARM and FOFEM) that simulate fire spread and effects based on fire spread physics developed from laboratory experiments [29,[33][34][35], and the empirical approaches that draw relationships from the analysis of existing fire severity patterns [36][37][38].In the past two decades, empirical relationships between fire severity pattern and its spatial controls have gained considerable attention [37][38][39][40], with scientific aims to identify the primary drivers and quantify their explanatory power on spatial heterogeneity of fire severity.
Since fire severity can be evaluated across a range of spatial scales, environmental variables in empirical models need to be consistently described at comparable spatial scales for maximum predictive power.Variables representing fuels and topography are considered to be bottom-up controls because their influences on fire patterns are pronounced largely at fine spatial scales [41].Those variables can be easily observed and/or resampled at various scales without losing their precision.For example, when fire severity is assessed in the field (site scale), the fuel properties and topography measured in situ are often used to build relationships with fire severity [39,40].When fire severity is quantified at coarser scales (e.g., fire patch level), vegetation variables derived from remotely sensed forest maps or forest inventory analysis databases, and topographic features developed from digital elevation models (DEMs) are often re-sampled to a coarser spatial resolution (upscaling) [42][43][44].Vegetation and topography have been widely recognized as the dominant controls on fire severity in many types of forest ecosystems [37,38,42,43,45].
Fire weather is considered a top-down control of wildfires as its influences on fire behavior are pronounced at broad scales.Although fire weather has been proven an influential driver of fire severity at the fire patch level [37,42,43,46,47], its influences and predictive power on spatial heterogeneity of fire severity remain poorly understood at fine scales.Current meteorological data are usually collected at a coarse spatial resolution, which is inconsistent with the fine scale at which site-level fire severity is assessed and predicted.The viewpoints that emphasize that fire weather is less important than vegetation or topography on regulating fire severity may be problematic as the low-resolution Forests 2018, 9, 130 3 of 26 fire weather may not represent the spatial variations of in situ meteorological conditions.Spatially interpolated data from weather stations can mitigate this issue, but in regions where the weather stations are scarce, this technique is inadequate.The lack of availability of credible and timely weather data limits the understanding of spatial controls on fire severity.
Recent advances in remote sensing techniques have provided a new opportunity to examine fine-scale fire weather effects on fire severity.Fire weather conditions are often characterized using the metric of fuel moisture content, which reflects dryness of dead fuels and water deficit of live biomass.Remotely sensed surface evapotranspiration (ET) products, which capture broad relationships between surface moisture availability and fuel dryness and vegetation drought-stress [44,48,49], have seldom been applied in modeling fire severity, even though they possess relatively high temporal and spatial resolutions.
Furthermore, prediction of fire severity patterns based on the empirical understanding of spatial controls has been well studied in North American boreal and western US forest landscapes [50][51][52], but there is still a lack of comprehensive analysis in Eurasian boreal forests, which are expected to become more fire prone with climate warming and drying [1,53,54].This study conducts a comprehensive analysis of within-patch fire severity variations in response to pre-fire vegetation, topography, and surface moisture availability at fine spatial scales in a Eurasian boreal forest landscape.The Great Xing'an boreal forest in northeastern China is an important forest ecosystem that stores 1.0-1.5 Pg C and provides 30% of the total timber yield in China [55].It is located near the southern frontier of Eurasian boreal forests where the fire regime is very sensitive to climate changes [43,54,56].Since April 2014, commercial logging has been completely forbidden in this region in an effort to restore and protect its valuable ecosystem services.Fire is the primary disturbance in this area, and there is an increasing urgency to understand the driving mechanisms of fire severity patterns to mitigate fire-induced ecological damage.In this paper, our objectives are three-fold.First, we investigate how fire severity varies across the landscape in response to environmental factors characterized at a consistent spatial resolution.Second, we verify the predictive power of remotely sensed pre-fire surface moisture conditions on determining fire severity patterns.Third, we develop an empirically based model to identify the distribution of potential high-severity burns.

Study Area
The boreal forest in the Great Xing'an Mountains of China is a fire-prone ecosystem that generally experiences frequent, moderate-to low-severity surface burns, mixed with infrequent high-severity crown fires.The climate is classified as mid-latitude continental cold-temperate type with short, warm, humid summers and long, cold, dry winters [28].The study area is mainly located in the Huzhong Forestry Bureau (Figure 1), which is situated in the central part of Great Xing'an region and represents a typical boreal forest landscape of northeastern China.It has a mean annual precipitation of approximately 460 mm that mostly occurs between July and September, and a mean annual temperature of approximately −4.7 • C. The topography is mountainous, with elevations ranging between 360 m and 1511 m above sea level.In contrast to the boreal forests dominated by evergreen coniferous tree species in North America and Europe, the forests in this area are dominated primarily by a deciduous coniferous tree species, Dahurian larch (Larix gmelinii (Rupr.)Rupr.), and mixed with some evergreen coniferous tree/shrub species including Korean spruce (Picea koraiensis Nakai), Scotch pine (Pinus sylvestris var.mongolica), and Siberian dwarf pine (Pinus pumila (Pall.)Regel), and a few deciduous broadleaf species of birch (Betula platyphylla) and aspen (Populus davidiana and Populus suaveolens Fisch.).The understory species are composed of evergreen shrubs (e.g., Ledum L. and Vaccinium vitis-idaea L.), deciduous shrubs (e.g., Betula fruticose Pall.and Rhododendron dauricum L.), and some herbaceous plants (e.g., Chamaenerion angustifolium (L.) and Carex appendiculata (Trautv.)Kukenth.)[28], whose distributions are influenced by the topographic and soil conditions [57].
Based on fire occurrence records published by the Chinese Forestry Science Data Center, there were 146 fires in the Huzhong Forestry Bureau's jurisdiction between 1991 and 2010, of which 111 were lightning-ignited [58].Most fires occurred in June, July, and August, which suggests that summer is the primary fire season in our study area.Similar fire occurrence tendencies for the entire Great Xing'an Mountains were reported from 1980 to 2005 in Fan et al. (2017) [54].Here we focused on 21 fires (Table 1) occurring in two fire years: 2000 (3 fires) and 2010 (18 fires).These two fire years exhibited the greatest burned areas of any other years in the past two decades, together accounting for about 82.5% of the total burned area over a 20-year period [28,43].Another important reason for selecting these fires is that they were all lighting-ignited in mid-to-late June and located within similar biotic and abiotic environments.In addition, field measurements of fire severity and forest regeneration in the area of these fires have been conducted by our research team since 2010 [28,59,60].  1) occurring in two fire years: 2000 (3 fires) and 2010 (18 fires).These two fire years exhibited the greatest burned areas of any other years in the past two decades, together accounting for about 82.5% of the total burned area over a 20-year period [28,43].Another important reason for selecting these fires is that they were all lighting-ignited in mid-to-late June and located within similar biotic and abiotic environments.In addition, field measurements of fire severity and forest regeneration in the area of these fires have been conducted by our research team since 2010 [28,59,60].

Remote Sensing Imagery Processing
We obtained four L-1 terrain-corrected Landsat TM and ETM+ images (path-row 121/24) with very good quality from 1999, 2000, 2007 and 2010 from the USGS website.To minimize spectral bias caused by phenology differences, the four Landsat images selected for study were all acquired in September.The raw digital number (DN) images of each spectral band were first calibrated into at-satellite radiance using the sensor-specific parameters cited in Chander et al. [61].A consistent radiometric response between multitemporal Landsat images is critically important for regional fire severity assessment over long time scales.Atmospherically-corrected Landsat surface reflectance products recently have been provided by the USGS, but to our knowledge, a further radiometric normalization can still be useful for consistent monitoring of forest changes.Thus, before using these Landsat images, we applied the 6S atmospheric correction method [62] and an absolute radiometric normalization approach, the iteratively reweighted multivariate alteration detection (IR-MAD) [63], to eliminate atmospheric effects and improve radiometric consistency between the Landsat time-series datasets.We selected a 2802 × 3483-pixel subset from each image (Figure 2), which covers all 21 fires, to carry out the normalization procedure in ENVI/IDL 4.7 (ITT Industries Inc., White Plains, NY, USA, 2009).We selected the 2007 image as the common reference image due to the least cloud coverage and used the bandwise regression parameters generated by IR-MAD to normalize the other three reflectance images.All the images were normalized to a consistent radiometric standard of the 2007 reference image.The dNBR is a well-known spectral index and has been proven to correlate well with fire severity in many types of forest ecosystems, including the study area of this research [28].We calculated the NBR and dNBR indices using Landsat bands 4 and 7 as proposed in Key and Benson [19]: Forests 2018, 9, 130 6 of 26 The two images acquired in 1999 and 2007 were used to calculate pre-fie NBR for the years of 2000 and 2010, respectively, while the images acquired in 2000 and 2010 were used to calculate post-fire NBR.
The 8-day moderate resolution imaging spectroradiometer (MODIS) ET product (MOD16A2) includes actual ET (AET), latent heat flux, potential ET (PET), and potential latent heat flux data.It is composed of daily canopy evaporation, plant transpiration, and soil evaporation and is calculated as the average value of cloud-free ET during an 8-day period [64].MOD16A2 products have been evaluated based on flux tower measurements in many ecosystems and exhibit agreement with the ground-measured ET in eastern Asian forests [64,65].In this study, we assumed that the MODIS-derived SMA index captures the broad relationship between remote sensing spectral signals and fuel moisture content.Although we did not develop an empirical relationship model to retrieve the actual pre-fire fuel moisture content, many recent publications have tested this hypothesis and indicated that MODIS data could be operationally integrated into fire danger systems [48,49,66].Currently, there are two versions of MOD16A2 products available.The latest (version 6) can provide 500 m ET observations, but it is not available for year 2000 as its collection began in 2001.Version 5 is available for both fire years, but its spatial resolution is coarser (1 km) (Figure A1).In addition, because these two versions used different model input datasets (e.g., land cover product, leaf area index) and algorithms [67], the final ET outputs differ between version 5 and version 6 (Figure A2).We obtained both version 5 and version 6 of MOD16A2 product for our study area and assessed the performance of both versions in our modeling.In our study, we mainly used the version 5 product for a consistent comparison of the models for 2000 and 2010.We also compared model performance and the ability of versions 5 and 6 in explaining spatial variation of fire severity for the year 2010.that MODIS data could be operationally integrated into fire danger systems [48,49,66].Currently, there are two versions of MOD16A2 products available.The latest (version 6) can provide 500 m ET observations, but it is not available for year 2000 as its collection began in 2001.Version 5 is available for both fire years, but its spatial resolution is coarser (1 km) (Figure A1).In addition, because these two versions used different model input datasets (e.g., land cover product, leaf area index) and algorithms [67], the final ET outputs differ between version 5 and version 6 (Figure A2).We obtained both version 5 and version 6 of MOD16A2 product for our study area and assessed the performance of both versions in our modeling.In our study, we mainly used the version 5 product for a consistent comparison of the models for 2000 and 2010.We also compared model performance and the ability of versions 5 and 6 in explaining spatial variation of fire severity for the year 2010.Based on occurrence dates of 21 fires (Table 1), we obtained five pre-fire MOD16A2 ET data (tile: H25V03, day of year from 129 to 168) for 2000 and six 1 km pre-fire MOD16A2 ET data (tile: H25V03, day of year from 129 to 176) for 2010 from the University of Montana's Numerical Terra-dynamic Simulation group.The six 500 m MOD16A2 ET data (tile: H25V03, day of year from 129 to 176) were obtained from the USGS Land Processes Distributed Active Archive Center.With the assistance of MOD16A2 Quality Control (QC) data, we only selected high-quality pixels (QC equal 0) from each 8-day ET dataset and used them to generate two integrated pre-fire ET datasets for 2000 and 2010.For a given burned pixel, the maximum of six observations with high quality were ranked by the time since fire occurrence date.We gave preference to the observation whose acquisition time was closest to the fire occurrence date to ensure the selection reflects the latest pre-fire surface ET conditions.Based on occurrence dates of 21 fires (Table 1), we obtained five pre-fire MOD16A2 ET data (tile: H25V03, day of year from 129 to 168) for 2000 and six 1 km pre-fire MOD16A2 ET data (tile: H25V03, day of year from 129 to 176) for 2010 from the University of Montana's Numerical Terra-dynamic Simulation group.The six 500 m MOD16A2 ET data (tile: H25V03, day of year from 129 to 176) were obtained from the USGS Land Processes Distributed Active Archive Center.With the assistance of MOD16A2 Quality Control (QC) data, we only selected high-quality pixels (QC equal 0) from each 8-day ET dataset and used them to generate two integrated pre-fire ET datasets for 2000 and 2010.For a given burned pixel, the maximum of six observations with high quality were ranked by the time since fire occurrence date.We gave preference to the observation whose acquisition time was closest to the fire occurrence date to ensure the selection reflects the latest pre-fire surface ET conditions.

Fire Severity Mapping
Burned pixels (30 m) were extracted based on a thresholding process for the forest disturbance index following the protocol of the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [68].The detailed description of this approach can be found in Fang et al. (2015) [43].Twenty-one fires burned about 22,650 ha in total; the three fire events in 2000 burned over 12,880 ha (about 251, 670 Landsat pixels), and the 18 fires in 2010 burned over 9760 ha (about 108, 450 Landsat pixels).These burned pixels exhibited significantly different spectral features compared to the unburned pixels in the Landsat images, as shown in Figure 2. Our previous study established a quadratic polynomial relationship between the dNBR and field measured composite burn index (CBI): which was confirmed to explain 84.6% variance in 74 CBI plots of 2010 fires and to produce accurate severity maps of the largest fire in 2000 (Kappa = 0.72) [28].The CBI values of 1.1 and 2.0 were applied as boundaries of three severity levels (i.e., low severity (0.1 ≤ CBI ≤ 1.1), moderate severity (1.1 < CBI ≤ 2), and high severity (2 < CBI ≤ 3)) as they were corresponding to about 10% and 80% of canopy mortality, respectively, in our ecosystem, where the live mature trees were very important for vegetation restoration [28].The two dNBR thresholds associated with CBI values of 1.1 and 2.0 were 0.484 and 1.099, respectively; these thresholds were applied to classify burned pixels of 2000 and 2010 into three severity levels.

Environmental Metrics
We used a suite of explanatory variables to describe various aspects of fuel, topography, and surface moisture conditions, and modeled their relationships with fire severity (Table 2).Using Landsat imagery, we developed two vegetation cover maps for two different years (1999 and 2007) and reconstructed the pre-fire vegetation conditions.Using a stratified decision tree classification method, we combined the conspicuous differences of phenology and spectral characteristics among vegetation types, and used Landsat surface reflectance, spectral indices, and components of Tasseled Cap transformation to produce vegetation cover maps at 30 m resolution [69].The whole area was classified into six vegetated categories consisting of deciduous coniferous forest (DCF, i.e., larch forest), deciduous broadleaf forest (DBF), evergreen coniferous forest (ECF), mixed forest (MF), grassland (GRS), and shrublands (SRB), as well as five non-vegetation categories (water bodies, bare rock, bare soil, urban land, and shade from mountains and clouds).Accuracy assessments were carried out in Fang et al. (2015) [43], which reported that the overall accuracy and the Kappa coefficient of the vegetated areas were 64% and 57% respectively.In detail, the mapping accuracy of DCF is 72%, DBF is 40%, ECF is 81%, MF is 49%, GRS is 82%, and SRB is 70%, indicating suitability for subsequent analysis.
We used a 30 m ASTER global digital elevation model (GDEM) product published by the National Aeronautics and Space Administration (NASA) to derive a number of commonly used topographic indices, including elevation (ELV), slope (SLP), aspect, and topographic wetness index (TWI).Aspect was further converted to better represent the potential influence of solar radiation (PSR) on site moisture conditions, following the equation: where θ is the aspect in degrees [70].The higher PSR values represent higher potential solar insolation.
The TWI was designed based on the assumption that water movement is controlled by the topography of the slopes [71].TWI is defined as: where α is the local upslope area draining through a certain point per contour length and tanβ is the slope angle at the point [72].TWI is often used to describe the spatial distribution of soil moisture and surface saturation, and has shown good correlation with soil moisture and depth to ground water in Eurasian boreal forests [72,73].TWI values typically range from 3 to 30, where a higher value indicates higher soil moisture potential.Fuel moisture content is a good surrogate for representing the indirect impacts of weather on fire behavior [74] because it influences many fire processes, such as ignition, combustion, and smoldering [75].Fuel moisture is also controlled by weather conditions such as precipitation and ET and is closely associated with soil moisture [48,75].Here, we calculated surface moisture availability (SMA) as the indicator of live fuel moisture content using the equation: where AET is actual ET and PET is the potential ET [76].The SMA is a critical parameter in governing the partition between sensible and latent heat flux at the surface, which is the key determinant of soil moisture, and also determines the water stress and water content of live plants [49,76].In general, when AET is equal to PET, the surface will reach moisture-saturated conditions and the SMA will be equal to 1.However, when the SMA is below certain threshold values, the vegetation may begin to suffer from drought stress and the soil may approach a limiting dryness [76].

Spatial Data Processing
To mitigate the issue of inconsistent spatial resolutions between various spatial datasets, we conservatively selected a 240 m spatial resolution as the overall standard.This resolution (240 m) is exactly eightfold the resolution (30 m) of the Landsat-derived vegetation maps and the ASTER GDEM data.This resolution is particularly suitable for aggregating the fine-resolution categorical map to the coarse-resolution continuous imagery.In addition, it is close to 250 m, which is the finest spatial resolution of MODIS products, and many remote sensing studies use 250 m to conduct downscaling for 1 km MODIS ET products [77].To downscale the SMA data, we used a nearest neighbor interpolation approach to resample the data from 1 km to 240 m.
Three different spatial aggregation methods were employed to upscale the vegetation, topographic, and fire severity data.Based on the Landsat-derived vegetation cover maps, we calculated the fractional coverage of each vegetation type by calculating the proportion of occurrence within an 8 × 8 pixel window (8-PW).Thus, the 30 m categorical vegetation maps were converted into 240 m continuous vegetation coverage images.For the numerical ASTER GDEM data, we calculated an average value for each 8-PW and subsequently used the aggregated DEM data to calculate the abovementioned topographic metrics.
To generate a fire severity value that can be consistently interpreted at different spatial scales, we used an area-weighted average method to aggregate the 30 m categorical fire severity map to 240 m.We first analyzed the proportions of unburned, low, moderate, and high severity pixels within each 8-PW; pixel proportions in the 8-PW serve as the area-weight values in the aggregation function.For instance, if 32 pixels within a given 8-PW (64 pixels) were classified as high severity, the weight value of high severity would be 50% (i.e., 32/64).If no pixel was classified as high severity in the 8-PW, the weight value of high severity would be 0.Then, a set of fire severity rating integers from 0 to 3 was used to indicate the fire severity gradient from unburned to high severity, respectively: where ω is the weight value of four severity gradients.The area-weighted average value of fire severity ranged from 0 to 3, with higher values representing higher fire severity.The area-weighted averaging process considers the relative importance of each severity level and provides a more balanced interpretation of the data.

Statistical Modeling
A boosted regression tree (BRT) model was applied to explain the relationships between remote sensed fire severity and selected environmental variables.The BRT model is a useful machine learning approach that uses recursive binary splits and a boosting technique to combine a large number of sequential trees to improve the fit and predictive performance [78,79].It has the capacity to handle complex relationships among numerical and categorical variables, quantify interactions between explanatory variables and overcome inaccuracies associated with regression and classification methods [78].We applied the "gbm" package (version 2.1.3)in R 3.4.1 (R Development Core Team 2017, Boston, MA, USA) to run the BRT analysis.We ran two sets of BRT models for the two different fire years, 2000 and 2010.We used the same parameter settings and amounts of training data to ensure that the outcomes of the two models were comparable.The random subsampling and bagging procedures in a BRT model may introduce stochasticity into the model outcomes.To mitigate the stochastic errors and create stable model outputs, we carried out 50 BRT modeling trials independently and calculated an average as the final result.
Before running the BRT models, we carried out a data exploration procedure to assess whether there was any conspicuous spatial autocorrelation or collinearity.Using functions in the "ape" package [80], we calculated the Moran's I index to examine spatial autocorrelation in fire severity and explanatory variables.Moran's I is similar to a correlation coefficient and represents the similarity of an observation to its nearby observations [81].It ranges between −1 and 1.Higher positive values indicate greater similarity, which can be interpreted as a spatially clustered distribution, while lower negative values indicate stronger dissimilarity (i.e., a more dispersed spatial pattern), and the zero value indicates a random spatial pattern (i.e., perfect randomness).We set a 300 m sample spacing distance (>1 pixel) when computing the Moran's I index.Although p-values were less than the 0.05 level, the Moran's I values of both response and explanatory variables were small (mostly less than 0.20, see Table A1), suggesting weak spatial autocorrelations among sampling pixels.Similar or higher Moran's I thresholds for determining suitable sampling spaces have been used in wildfire-related literature [36,45,50].We also calculated the pairwise Pearson correlation coefficient (r) using the "Hmisc" package to evaluate potential collinearity among predictor variables.All variables selected for modeling had low pairwise Pearson correlations (|r| < 0.60), suggesting relatively low levels of collinearity.We used a subset of 200 random sampling pixels to build the BRT model, while the remaining 100 pixels were used for validation.
The optimization of a BRT model is jointly controlled by the learning rate, tree complexity, bagging fraction, and number of trees.In this study, we set these parameters at 0.01, 5 and 0.5, respectively, following the recommended model inputs of Elith et al. (2008) [79].The number of trees in each BRT trial was automatically selected based on a 5-fold cross-validation procedure to avoid over-fitting problem.The relative importance of each explanatory variable to fire severity was measured by averaging the frequency it was selected for splitting among all trees, and the importance was weighted by the squared improvement to the model as a result of each split [79].The relative importance values of the predictors were scaled as a percentage, the sum of which equals 100%.Higher importance values represented stronger impacts on controlling fire severity.In addition, we explored the dependency relationships between several important variables and fire severity by plotting the effect of a specific explanatory variable on the response variable after averaging the effects of remaining explanatory variables in the same model.
The coefficient of determination (R 2 ) reported by the BRT model was used to evaluate how well the model fits the training data.We used the 100 samples that were independent of the training data to assess the predictive power of our BRT models.We examined the goodness of fit between simulated fire severity values and observed fire severity values by applying a linear regression model and calculating the squared multiple correlation coefficients (denoted r 2 ).To provide a standalone evaluation of model performance, we used the error matrix method to evaluate the predicted fire severity maps.We selected an optimal simulation image for each fire year based on the r 2 value and classified the burned pixels into low (fire severity ≤ 1), moderate (1 < fire severity ≤ 2), and high (2 < fire severity ≤ 3) severity levels.The two aggregated fire severity images were also converted into categorical severity maps using these same thresholds and were then used as reference maps.We used all burned pixels from each fire severity level to ensure unbiased evaluation for both fire years.

Evaluation of Model Performance
We found that the two sets of BRT models fit the training data very well, as they explained 83.0% (R 2 max = 85.0%,SD = 0.008) of the variation in fire severity in 2000 (Figure 3a) and 73.8% (R 2 max = 81.3%,SD = 0.035) of the variation in fire severity in 2010 (Figure 3b).Similarly, the 50 validation models of the two years showed goodness of fits of 0.795 (r 2 max = 0.820, SD = 0.040) (Figure 3c) and 0.618 (r 2 max = 0.656, SD = 0.012) (Figure 3d), respectively, suggesting that the predicted fire severity results were well correlated with the reference values.
The confusion matrices indicated that the optimal BRT model in 2000 achieved higher prediction accuracy than the optimal model in 2010 (Table 3).The predicted severity map of 2000 produced fewer commission errors in moderate (44.3% vs. 71.6%)and high (11.2% vs. 24.0%)severity pixels than 2010, but it generated a high commission error for low severity (85.8%) pixels, indicating more pixels were erroneously predicted as low severity than were actually observed in the severity results.The predicted severity map of 2010 generated high commission error (71.6%) for moderate severity level and high omission error (64.6%) for low severity level, indicating overestimated moderate severity area and underestimated low and high severity area (Figure A3).Overall, we found the prediction map of 2000 to represent very good agreement with the observed severity map, while the prediction map of 2010 underestimated the severity of most of the 18 fires (Figure 4).

Relative Importance of Environmental Variables
The relative importance of individual environmental variables varied substantially in the different fire years (Figure 5).For the fires in 2000, the six most important predictors of fire severity in decreasing order were DBF, SRB, MF, TWI, SLP and ELV.These predictors contributed to a total of 76.8% of the relative importance.Three vegetation variables together contributed over 55% relative importance, indicating the spatial pattern of fire severity in 2000 was largely driven by the distribution of these three vegetation types.For the 2010 burns, the six most important predictors were ECF, SMA, SLP, PSR, ELV and TWI, which together contributed to 75.1% of the relative importance.The ECF and SMA variables independently contributed as much as 23.3% and 12.8% relative importance respectively.The 500 m version 6 SMA was found to have a slightly lower (11.0%)

Relative Importance of Environmental Variables
The relative importance of individual environmental variables varied substantially in the different fire years (Figure 5).For the fires in 2000, the six most important predictors of fire severity in decreasing order were DBF, SRB, MF, TWI, SLP and ELV.These predictors contributed to a total of 76.8% of the relative importance.Three vegetation variables together contributed over 55% relative importance, indicating the spatial pattern of fire severity in 2000 was largely driven by the distribution of these three vegetation types.For the 2010 burns, the six most important predictors were ECF, SMA, SLP, PSR, ELV and TWI, which together contributed to 75.1% of the relative importance.The ECF and SMA variables independently contributed as much as 23.3% and 12.8% relative importance respectively.The 500 m version 6 SMA was found to have a slightly lower (11.0%)relative importance than the 1 km version 5 SMA product when modeling spatial variability of fire severity in 2010 (Figure A4).An overall ranking of 11 variables was calculated based on a weighted average approach (see y-axis of Figure 5), which identified four vegetation variables (ECF, DBF, SRB, and MF) as the primary controls of fire severity, followed by SMA and topographic variables.relative importance than the 1 km version 5 SMA product when modeling spatial variability of fire severity in 2010 (Figure A4).An overall ranking of 11 variables was calculated based on a weighted average approach (see y-axis of Figure 5), which identified four vegetation variables (ECF, DBF, SRB, and MF) as the primary controls of fire severity, followed by SMA and topographic variables.relative importance than the 1 km version 5 SMA product when modeling spatial variability of fire severity in 2010 (Figure A4).An overall ranking of 11 variables was calculated based on a weighted average approach (see y-axis of Figure 5), which identified four vegetation variables (ECF, DBF, SRB, and MF) as the primary controls of fire severity, followed by SMA and topographic variables.When the 11 variables were grouped into three broad control types of vegetation, topography, or SMA, the relative importance of these three control types became quite similar between the two fire years.Vegetation consistently played a strong role in determining fire severity in the two fire years because it contributed about 40-80% of total relative importance (Figure 6a).Together with 20-50% of the total relative importance contributed by topographic variables, the results revealed that potential fire severity in our study area is mostly controlled by vegetation and topography.At the same time, we found SMA has similar relative importance to the maximum importance value of topographic variables, indicating that SMA also has considerable predictive power over fire severity (Figure 6b).by multiplying two mean R 2 (i.e., 0.830 for RIMPs of 2000, 0.738 for RIMPs 2010) with two mean RIMP values (i.e., RIMP2000 and RIMP2010) of each variable.The ECF has the highest overall relative importance while the GRS has the lowest value.See Table 2 for definition of variable abbreviations.
When the 11 variables were grouped into three broad control types of vegetation, topography, or SMA, the relative importance of these three control types became quite similar between the two fire years.Vegetation consistently played a strong role in determining fire severity in the two fire years because it contributed about 40-80% of total relative importance (Figure 6a).Together with 20-50% of the total relative importance contributed by topographic variables, the results revealed that potential fire severity in our study area is mostly controlled by vegetation and topography.At the same time, we found SMA has similar relative importance to the maximum importance value of topographic variables, indicating that SMA also has considerable predictive power over fire severity (Figure 6b).

Relationships between Environmental Variables and Fire Severity
Partial dependence plots were helpful for visualizing the response of fire severity to explanatory variables.We found that the relationships between fire severity and environmental variables were mostly nonlinear and varied little in different years (Figure 7).In general, high severity fires were more likely to occur in dense evergreen forests (Figure 7a), especially on well-drained, gentle slopes at high altitudes (Figure 7f-h).On the other hand, the partial dependence of fire severity generally decreased with increasing DBF and SRB (Figure 7b and c), suggesting that the fire-resistant traits of these two vegetation types may mitigate the adverse effects of severe burning.We found that increased coverage of mixed forests generally had a negative impact on fire severity (Figure 7d), especially in the year 2000, in which MF was recognized as an important variable.Fire severity in deciduous coniferous forests was significantly lower than in evergreen coniferous forests (Figure 7j vs. 7a) but higher than in the other forest types.DCF did not exhibit a strong relationship with fire severity in 2010, but did exhibit a positive relationship with fire severity in 2000.The GRS variable was found to have little relative importance to fire severity and represented a weak negative impact on fire severity in 2010 (Figure 7k).The south facing slopes and areas with high TWI values typically experienced moderate severity fires (Figure 7h,i).Our results also indicated that higher probabilities

Relationships between Environmental Variables and Fire Severity
Partial dependence plots were helpful for visualizing the response of fire severity to explanatory variables.We found that the relationships between fire severity and environmental variables were mostly nonlinear and varied little in different years (Figure 7).In general, high severity fires were more likely to occur in dense evergreen forests (Figure 7a), especially on well-drained, gentle slopes at high altitudes (Figure 7f-h).On the other hand, the partial dependence of fire severity generally decreased with increasing DBF and SRB (Figure 7b,c), suggesting that the fire-resistant traits of these two vegetation types may mitigate the adverse effects of severe burning.We found that increased coverage of mixed forests generally had a negative impact on fire severity (Figure 7d), especially in the year 2000, in which MF was recognized as an important variable.Fire severity in deciduous coniferous forests was significantly lower than in evergreen coniferous forests (Figure 7j vs. Figure 7a) but higher than in the other forest types.DCF did not exhibit a strong relationship with fire severity in 2010, but did exhibit a positive relationship with fire severity in 2000.The GRS variable was found to have little relative importance to fire severity and represented a weak negative impact on fire severity in 2010 (Figure 7k).The south facing slopes and areas with high TWI values typically experienced moderate severity fires (Figure 7h,i).Our results also indicated that higher probabilities of moderateto high-severity fires were related to canopy or surface dryness (Figure 7e), especially when the 1 km SMA ranged from 0.13 to 0.30 in the summer.In contrast, we found a strong positive association between fire severity and 500 m SMA (Figure A5), as a result of inconsistent SMA values.
of moderate-to high-severity fires were related to canopy or surface dryness (Figure 7e), especially when the 1 km SMA ranged from 0.13 to 0.30 in the summer.In contrast, we found a strong positive association between fire severity and 500 m SMA (Figure A5), as a result of inconsistent SMA values.2.

Environmental Influences on Fire Severity
We found that the BRT method was effective in investigating the relationships between fire severity and environmental gradients, such as vegetation composition, terrain, and surface moisture status.Fire severity is a complex function of these environmental gradients and such relationships may vary in different years and locations.With thousands of 240 m sampling data from representative historical fires, we found that the spatial distribution of fire severity could be predicted with adequate precision in a Great Xing'an boreal forest landscape.In our exploration of the relative importance of these spatial controls on fire severity, we found that fuel conditions are the most influential predictor in determining the magnitude of fire severity.This finding is in accordance with our previous study [43], which was conducted at burned patch level for the entire Great Xing'an boreal forest, and studies in similar ecosystems, such as Canadian boreal forests [82,83] and subalpine forests [37,84,85].2.

Environmental Influences on Fire Severity
We found that the BRT method was effective in investigating the relationships between fire severity and environmental gradients, such as vegetation composition, terrain, and surface moisture status.Fire severity is a complex function of these environmental gradients and such relationships may vary in different years and locations.With thousands of 240 m sampling data from representative historical fires, we found that the spatial distribution of fire severity could be predicted with adequate precision in a Great Xing'an boreal forest landscape.In our exploration of the relative importance of these spatial controls on fire severity, we found that fuel conditions are the most influential predictor in determining the magnitude of fire severity.This finding is in accordance with our previous study [43], which was conducted at burned patch level for the entire Great Xing'an boreal forest, and studies in similar ecosystems, such as Canadian boreal forests [82,83] and subalpine forests [37,84,85].
Coniferous forests/shrubs experience higher severity fires than broadleaf forests and shrublands in our ecosystem.Such result is generally consistent with observations in North America boreal forests where deciduous tree dominated stands are found to be fire break and reduce landscape flammability owing to higher foliage moisture and less surface fuels [86,87].Larix gmelinii is the dominant coniferous tree species in the boreal forests of Northeastern China (Figure 8a).Unlike the dominant evergreen coniferous species (e.g., black spruce) in Northern American boreal forests or subalpine forests (e.g., spruce-fir, ponderosa pine) in Western US, larch is a deciduous coniferous species; it is considered fire-tolerant because its self-pruning trait can reduce ladder fuels, and the open crown can reduce the bulk density and connectivity of canopy fuels [88].However, we found that larch forests usually experience high mortality in flat areas and valley bottoms due to heat-induced root damage (Figure 8c).To adapt to shallow soil and permafrost, Larix gmelinii develops lateral roots in the thick, flammable moss and organic soil layers to improve adhesion and nutrient supply, but this also means it can be easily injured or killed by surface fires [43].Larch forests are also characterized by abundant understory plants due to their open-canopy environment, which can contribute higher fire severity than the forests without an abundant understory component.Another important coniferous species in this region is Pinus pumila (Figure 8b), which is evergreen and usually mixed with Larix gmelinii in open forests at altitudes of 800-1200 m or grows densely on rocky ridges at altitudes higher than 1200 m [89].Pinus pumila is highly flammable because it contains abundant volatile organic compounds in its needles, twigs, and seeds [90].Furthermore, windy and dry conditions can accelerate the spreading of fires on the ridges of mountains.Therefore, as demonstrated in our analysis, the increases in Pinus pumila coverage may considerably increase the probability of high severity fires (Figure 8d).Such finding is in consistent with Estes et al. (2017) [38], who reported that shrub species that are favored by fires can generate higher fire severity than mixed hardwood/coniferous forests and hardwood forests in northern California.Coniferous forests/shrubs experience higher severity fires than broadleaf forests and shrublands in our ecosystem.Such result is generally consistent with observations in North America boreal forests where deciduous tree dominated stands are found to be fire break and reduce landscape flammability owing to higher foliage moisture and less surface fuels [86,87].Larix gmelinii is the dominant coniferous tree species in the boreal forests of Northeastern China (Figure 8a).Unlike the dominant evergreen coniferous species (e.g., black spruce) in Northern American boreal forests or subalpine forests (e.g., spruce-fir, ponderosa pine) in Western US, larch is a deciduous coniferous species; it is considered fire-tolerant because its self-pruning trait can reduce ladder fuels, and the open crown can reduce the bulk density and connectivity of canopy fuels [88].However, we found that larch forests usually experience high mortality in flat areas and valley bottoms due to heatinduced root damage (Figure 8c).To adapt to shallow soil and permafrost, Larix gmelinii develops lateral roots in the thick, flammable moss and organic soil layers to improve adhesion and nutrient supply, but this also means it can be easily injured or killed by surface fires [43].Larch forests are also characterized by abundant understory plants due to their open-canopy environment, which can contribute higher fire severity than the forests without an abundant understory component.Another important coniferous species in this region is Pinus pumila (Figure 8b), which is evergreen and usually mixed with Larix gmelinii in open forests at altitudes of 800-1200 m or grows densely on rocky ridges at altitudes higher than 1200 m [89].Pinus pumila is highly flammable because it contains abundant volatile organic compounds in its needles, twigs, and seeds [90].Furthermore, windy and dry conditions can accelerate the spreading of fires on the ridges of mountains.Therefore, as demonstrated in our analysis, the increases in Pinus pumila coverage may considerably increase the probability of high severity fires (Figure 8d).Such finding is in consistent with Estes et al. (2017) [38], who reported that shrub species that are favored by fires can generate higher fire severity than mixed hardwood/coniferous forests and hardwood forests in northern California.Our results indicated that topographic factors had a considerable influence on fire severity.It is well known that topography can influence fire behavior by impacting fuel moisture, local wind patterns, fire spread direction, and vegetation composition [32], but quantitative demonstrations of these relationships are still needed for optimal mitigation of adverse fire effects.Our results showed that slope and elevation are the two most important topographic variables, followed by TWI and PSR.This suggests that the primary pathway by which topography regulates fire severity in this Siberian boreal ecosystem is by governing fuel moisture and by strongly interacting with vegetation conditions.For example, we found that severe fires were more likely to occur in high altitude regions.A similar pattern was found in dry ponderosa pine forests of the western US and boreal forests in Europe [42,51,91].Surface fuels on the upper slopes could dry quickly due to efficient drainage and greater degrees of solar exposure, which may increase the flammability of fuels and facilitate more severe burns.The preheating effects of upslope fires on the adjacent fuels can also increase fire severity [32].In addition, it cannot be ignored that the ECF are generally densely distributed at high altitudes in this region.In general, north-facing slopes possess higher fuel moisture and lower surface temperature than south-facing slopes.However, in our ecosystem, we found that the north-facing slopes burned more severely due to higher biomass coverage [60].Similar findings were also reported in studies conducted in boreal and subalpine ecosystems [42,43,84].
Fuel moisture interacts with many complex ecological and physical processes, making it important yet difficult to represent in a modeling framework used to study its spatiotemporal dynamics and influences on fire regime [75].Topographic gradients could provide a partial explanation for the spatial variation of fuel moisture, especially in light of our finding that TWI is inversely associated with fire severity, as expected.Previous studies proved that TWI is closely associated with soil moisture in European boreal forests [72,92].Thus, we believe the drainage condition largely determines the moisture gradient of a forest stand and further influences fuel moisture dynamics.Dead fuel moisture dynamics are driven by three mechanisms-capillary forces, infiltration, and diffusion-among which, infiltration and diffusion are the primary driving mechanisms and are both influenced by moisture gradients [75].Although live fuel moisture is driven by different mechanisms than dead fuel moisture dynamics, soil water dynamics are still an important part of those mechanisms and can directly influence plant transpiration.
The MODIS-derived 1 km SMA exerted considerable influences on model performance, and its relationship with fire severity was negative, as expected.This finding aligns with the work of van Mantgem et al. (2013) [18], which suggests that a pre-fire water deficit can increase fire severity (tree mortality) because the drought-stressed trees are vulnerable to fire-induced injury.Similarly, Xiao and Zhuang [93] found that drought directly affected fire activity in Canadian and Alaskan boreal forests by enhancing fuel flammability and increasing ignitions.However, it should be noted that the pre-fire SMA applied in this study can only represent the short-term temporal variability of the surface moisture conditions, which may not necessarily reflect the long-term effects of drought stress on fire severity.It has been reported that plant communities within a forest stand, especially understory vegetation layers, may be influenced by the long-term drought stress that is regulated by topographic and climatic factors [94].
Although the ranking of relative importance of vegetation, topography, and SMA was similar between the two fire years (Figure 6), the relative importance of individual explanatory variables differed between the two fire years (Figure 5).Such differences may be attributable to their different pre-fire vegetation composition, structure, and disturbance history.By comparing the pre-fire vegetation compositions (Figure 9), we found that fires in 2000 had lower proportions of DBF and MF but higher proportions of DCF than the 2010 fires.Fires in 2000 burned more areas in the Huzhong Natural Reserve where it is dominated by mature (>100 years) larch forest (DCF) due to the strictly enforced cutting ban, and most DBF, MF, and shrubs are located in recently burned/disturbed areas that carry significantly less fuels.Consequently, the proportion of DBF, MF, and shrubs were considered more important than proportion of DCF in modeling severity of 2000 fires.In contrast, fires in 2010 burned more areas in the Huzhong Forestry Bureau jurisdiction that had been disturbed by clear cutting since 1950s, and as recent as 2000s, leading to greater abundance of young stands irrespective of forest type [95].This could partly explain why fire severity of 2010 in the areas with high proportions of DCF was similar to the areas with high proportion of DF or MF (Figure 7j vs. Figure 7b,d).In contrast, because the highly flammable evergreen coniferous shrub species Pinus pumila is not an economically viable species to cut, its fuel loading is generally higher than young stands of other forest types.Consequently, the fire severity in areas with high proportion of ECF was high (Figure 7a), and ECF was considered a more important vegetation variable in driving overall fire severity variability of 2010 fires.
irrespective of forest type [95].This could partly explain why fire severity of 2010 in the areas with high proportions of DCF was similar to the areas with high proportion of DF or MF (Figure 7j vs. 7b and 7d).In contrast, because the highly flammable evergreen coniferous shrub species Pinus pumila is not an economically viable species to cut, its fuel loading is generally higher than young stands of other forest types.Consequently, the fire severity in areas with high proportion of ECF was high (Figure 7a), and ECF was considered a more important vegetation variable in driving overall fire severity variability of 2010 fires.Variable abbreviations were described in Table 2.

Prediction of Fire Severity
Our results revealed that landscape-level fire severity is primarily determined by fuel and terrain features, suggesting that fire severity is predictable at this scale.The prediction accuracy is determined by the quality of the relevant predictors or proxies.Fuel conditions are usually complicated and difficult to characterize, not only because of the endogenous variety of plant communities throughout the landscape but also because of exogenous factors, such as disturbance regime, climate, and anthropogenic activities [75].Although various remote sensing methods were proposed to improve fuel mapping, there are still challenges in accurately quantifying the critical fuel parameters that can regulate fire effects, such as fuel loading and canopy bulk density [96,97].Remotely sensed spectral information is usually applied as ancillary data or proxies when modeling fuel parameters, as it is strongly correlated with many biophysical vegetation parameters, such as biomass, leaf area index, and productivity [26,98,99].
However, remotely sensed imagery cannot detect surface fuels obscured by the forest canopy and insufficiently distinguishes fine fuels from dead biomass pools due to the inconsistency between particle size and the spatial resolution of the image [100].Compared with the spatial complexity and high variability of fuel parameters at the landscape scale, vegetation cover type is relatively identifiable due to the unique spectral features of plant communities.The Landsat-derived vegetation cover data can reflect the general variability of vegetation composition but may not necessarily depict parameters related to fuel type or fuel particle size.Our results indicated that vegetation coverage could reliably explain the variability of fire severity at a 240 m spatial scale.The classification accuracy of vegetation mapping could also affect the predictability of fire severity models.The spatial aggregation process may confound the accuracy of the 240 m vegetation coverage data, thus increasing the uncertainty of the modeling results.
The large-scale terrain features of wildlands are usually invariant over long periods.In addition, terrain features can be accurately characterized using various traditional or modern survey technologies.Together with their close relationship with fire behavior, topographic variables are commonly used as predictors of fire severity.Digital elevation models provide essential topographic information from which aspect, slope, and other terrain features can be derived.However, it should be noted that some topographic conditions are inherently scale-dependent.Thus, the upscaling Variable abbreviations were described in Table 2.

Prediction of Fire Severity
Our results revealed that landscape-level fire severity is primarily determined by fuel and terrain features, suggesting that fire severity is predictable at this scale.The prediction accuracy is determined by the quality of the relevant predictors or proxies.Fuel conditions are usually complicated and difficult to characterize, not only because of the endogenous variety of plant communities throughout the landscape but also because of exogenous factors, such as disturbance regime, climate, and anthropogenic activities [75].Although various remote sensing methods were proposed to improve fuel mapping, there are still challenges in accurately quantifying the critical fuel parameters that can regulate fire effects, such as fuel loading and canopy bulk density [96,97].Remotely sensed spectral information is usually applied as ancillary data or proxies when modeling fuel parameters, as it is strongly correlated with many biophysical vegetation parameters, such as biomass, leaf area index, and productivity [26,98,99].
However, remotely sensed imagery cannot detect surface fuels obscured by the forest canopy and insufficiently distinguishes fine fuels from dead biomass pools due to the inconsistency between particle size and the spatial resolution of the image [100].Compared with the spatial complexity and high variability of fuel parameters at the landscape scale, vegetation cover type is relatively identifiable due to the unique spectral features of plant communities.The Landsat-derived vegetation cover data can reflect the general variability of vegetation composition but may not necessarily depict parameters related to fuel type or fuel particle size.Our results indicated that vegetation coverage could reliably explain the variability of fire severity at a 240 m spatial scale.The classification accuracy of vegetation mapping could also affect the predictability of fire severity models.The spatial aggregation process may confound the accuracy of the 240 m vegetation coverage data, thus increasing the uncertainty of the modeling results.
The large-scale terrain features of wildlands are usually invariant over long periods.In addition, terrain features can be accurately characterized using various traditional or modern survey technologies.Together with their close relationship with fire behavior, topographic variables are commonly used as predictors of fire severity.Digital elevation models provide essential topographic information from which aspect, slope, and other terrain features can be derived.However, it should be noted that some topographic conditions are inherently scale-dependent.Thus, the upscaling process used for the ASTER GDEM data may filter out some detailed terrain features that would be reflected at the 30 m spatial resolution.Because we focused on the 240 m spatial resolution, we did not examine the sensitivity of the relationships between fire severity and topographic variables to the scaling process.Many studies have reported scale dependency in the relationships between topographic characteristics Forests 2018, 9, 130 18 of 26 and fire attributes, such as severity, frequency, burned area, and burn probability [101][102][103][104].Despite the good performance of our BRT models, we conclude that 240 m may not be the optimal spatial resolution for predicting fire severity.Enlarging spatial scales (e.g., from 30 to 240 m) can be beneficial in refining the relationships between fire severity and environmental gradients, but it may decrease the visual effect of the prediction maps due to the coarse spatial resolution.Therefore, we suggest that future applications should weigh the model performance against the practical applicability of the prediction maps.

Limitation and Uncertainty
Despite efforts to improve the prediction of fire severity by incorporating sound explanatory variables, some knowledge gaps should be noted as they may influence the interpretation of the modeling results.First, each MOD16A2 pixel contains the best possible daily ET estimation during the 8-day period, which was selected based on the imaging conditions and observation coverage.The arbitrary application of SMA as a proxy for fuel moisture may increase the alternative quantification of land surface status for modeling fire severity, but its real relationship with actual fuel moisture still needs to be validated in our ecosystem.Due to the lack of daily fire progression maps, and to very sparse weather station coverage in our study area, we cannot address how day-to-day weather impacts fire severity in this study.With the advantage of high-frequency MODIS observations, many recent studies have begun to incorporate spatial interpolation approaches using MODIS data to characterize daily progressions of large fires [105,106].We believe such efforts can greatly improve the understanding of spatial controls on fire severity.For example, based on Landsat-derived fire progression maps and fire weather observations at 4-km spatial resolution, Birch et al. (2015) [37] investigated the influences of vegetation, topography, and daily fire weather on severity patterns of wildfires in the Western US and reported that vegetation cover had the greatest influence on fire severity; this is quite similar to our findings in this study, as well as in our previous patch-scale analysis [43].They also acknowledged that the coarse weather conditions may not fully reflect the influences of microscale meteorological conditions on severity patterns.The inconsistent temporal and spatial resolution among daily weather observations, vegetation, and topography can obscure the real effects of weather on fire severity.
Fire severity is a result of accumulated fire effects on forest ecosystems because the thick and moist organic layers in boreal regions can prolong the fire duration.Although we tried to balance the spatial resolution among explanatory variables, it is somewhat challenging to reflect environmental gradients and fire activities at both fine spatial and temporal scales using the data sources available to us.At the same time, because MODIS ET product updates led to considerable changes in SMA values, we found that 500 m SMA has different relative importance and influences on fire severity.Our intentions were not to arbitrarily justify which kind of ET product is more reliable for fire severity prediction, especially without sufficient validation of 500 m ET products with site-based flux observation, but any improvement in spatial resolution is valuable and further efforts are encouraged to verify the suitability of these products for specific ecosystems.
Although our study area, the Huzhong Forestry Bureau jurisdiction, is a representative forest landscape of Great Xing'an Mountains and shares a similar fire regime as other nearby areas [54,56], we could not conclude that fire severity patterns of the entire region are following the same causal mechanism at finer scales.The purpose of this study is not to establish a global prediction model for fire severity that can be generalized for all Chinese boreal ecosystems.There were fires occurring in meadow and wetland ecosystems, as well as in forests dominated by broadleaf trees in the southern part of the Great Xing'an region.We believe our results may not be suitable for predicting the severity of those kinds of wildfires.In addition, sampling data were extracted from fires occurring in summer; although the vegetation coverage and topographic variables adopted in this study are insensitive to intra-annual variability, the relationship between SMA and fire severity may change in different seasons and should be further investigated.

Conclusions
Although the parallel comparison of the two models did not show strictly consistent modeling relationships, the models generally demonstrated that fire severity was strongly controlled by the coverage of certain vegetation types that have high flammability or fire resistance.The topographic conditions can help determine the distribution of flammable plant types and communities.Topography can also directly influence fuel moisture and create firebreaks through the drainage systems.Remotely sensed fuel moisture proxies (such as MODIS ET products) were also proven to play important roles in modeling fire severity.These findings reveal that fire severity is predictable at the landscape scale in our study area, and its prediction can be improved by incorporating spatial variables related to fire behavior.Our study provides an overview of the hotspot areas within the landscape where severe fires are most likely distributed.Such mapping capabilities allow managers to optimize fuel treatment strategies by considering the vegetation, topography, and spatial patterns of land surface moisture.The modeling framework employed in our study can readily incorporate new observations and simulated spatial datasets, promoting the more reliable prediction of fire severity in the future.
intra-annual variability, the relationship between SMA and fire severity may change in different seasons and should be further investigated.

Conclusions
Although the parallel comparison of the two models did not show strictly consistent modeling relationships, the models generally demonstrated that fire severity was strongly controlled by the coverage of certain vegetation types that have high flammability or fire resistance.The topographic conditions can help determine the distribution of flammable plant types and communities.Topography can also directly influence fuel moisture and create firebreaks through the drainage systems.Remotely sensed fuel moisture proxies (such as MODIS ET products) were also proven to play important roles in modeling fire severity.These findings reveal that fire severity is predictable at the landscape scale in our study area, and its prediction can be improved by incorporating spatial variables related to fire behavior.Our study provides an overview of the hotspot areas within the landscape where severe fires are most likely distributed.Such mapping capabilities can allow managers to optimize fuel treatment strategies by considering the vegetation, topography, and spatial patterns of land surface moisture.The modeling framework employed in our study can readily incorporate new observations and simulated spatial datasets, promoting the more reliable prediction of fire severity in the future.2.

Figure 1 .
Figure 1.Location of study area (a) showing severity of wildfires occurring in 2000 (blue perimeters) and 2010 (pink perimeters).Most of the study area is located in the Huzhong Forestry Bureau (b) in the middle of Great Xing'an boreal forests (green patch in b and c) which administratively belongs to Heilongjiang province in northeastern China (c).One fire was located in E'lunchun County (a), which belongs to the Inner Mongolian part of Great Xing'an boreal forests.Forests within the Huzhong Natural Reserve are primarily natural forests because of a strictly enforced ban on commercial and salvage logging within the reserve since 1958, while forests outside the natural reserve experienced severe cutting since the 1950s.

Figure 1 .
Figure 1.Location of study area (a) showing severity of wildfires occurring in 2000 (blue perimeters) and 2010 (pink perimeters).Most of the study area is located in the Huzhong Forestry Bureau (b) in the middle of Great Xing'an boreal forests (green patch in b and c) which administratively belongs to Heilongjiang province in northeastern China (c).One fire was located in E'lunchun County (a), which belongs to the Inner Mongolian part of Great Xing'an boreal forests.Forests within the Huzhong Natural Reserve are primarily natural forests because of a strictly enforced ban on commercial and salvage logging within the reserve since 1958, while forests outside the natural reserve experienced severe cutting since the 1950s.

Figure 2 .
Figure 2. Pre-fire (a) and post-fire (b) false color Landsat TM images (R-TM5, G-TM4, and B-TM3) of the study area.The light pink patches in (b) indicate fires occurred in 2000, while dark red patches indicate fires occurred in 2010.

Figure 2 .
Figure 2. Pre-fire (a) and post-fire (b) false color Landsat TM images (R-TM5, G-TM4, and B-TM3) of the study area.The light pink patches in (b) indicate fires occurred in 2000, while dark red patches indicate fires occurred in 2010.

Forests 2018, 9 , 26 Figure 3 .
Figure 3. Two examples of linear relationships used for validating model performance of 2000 (a,c) and 2010 (b,d) based on 200 training samples (a,b) and 100 verification samples (c,d).Coefficients of determination (R 2 for training samples), and the squared multiple correlation coefficients (r 2 for verification samples) are also plotted.Blue solid lines show predicted linear regression fit.Black dashed lines represent 1:1 line.

Figure 3 .
Figure 3. Two examples of linear relationships used for validating model performance of 2000 (a,c) and 2010 (b,d) based on 200 training samples (a,b) and 100 verification samples (c,d).Coefficients of determination (R 2 for training samples), and the squared multiple correlation coefficients (r 2 for verification samples) are also plotted.Blue solid lines show predicted linear regression fit.Black dashed lines represent 1:1 line.

Figure 4 .
Figure 4. Observed fire severity (a) aggregated from Landsat observations versus the modeled fire severity (b) for 21 fires based on boosted regression tree models.All fire severity images were plotted at 240 m spatial resolution.

Figure 5 .Figure 4 .
Figure 5. Relative importance proportion (RIMP) of explanatory variables for 50 boosted regression tree models (Mean + SD) of fire severity in 2000 and 2010.The y-axis (ECF to GRS) shows an overall rank of 11 variables in descending order according to a weighted average of RIMP, which is calculated

Figure 4 .
Figure 4. Observed fire severity (a) aggregated from Landsat observations versus the modeled fire severity (b) for 21 fires based on boosted regression tree models.All fire severity images were plotted at 240 m spatial resolution.

Figure 5 .Figure 5 .
Figure 5. Relative importance proportion (RIMP) of explanatory variables for 50 boosted regression tree models (Mean + SD) of fire severity in 2000 and 2010.The y-axis (ECF to GRS) shows an overall rank of 11 variables in descending order according to a weighted average of RIMP, which is calculated

Figure 6 .
Figure 6.Ternary plot (a) showing the total relative importance proportion (RIMP) of three groups of explanatory variables (i.e., Vegetation, Topography and Surface Moisture Availability (SMA)) in controlling fire severity of 2000 (blue) and 2010 (red).In considering the difference of variable numbers among three groups, a variable which has the maximum RIMP (MaxRIMP) of each group is selected and compared (b).Axis in ternary plot (b) representing relative proportions of MaxRIMP for three types of variables.

Figure 6 .
Figure 6.Ternary plot (a) showing the total relative importance proportion (RIMP) of three groups of explanatory variables (i.e., Vegetation, Topography and Surface Moisture Availability (SMA)) in controlling fire severity of 2000 (blue) and 2010 (red).In considering the difference of variable numbers among three groups, a variable which has the maximum RIMP (MaxRIMP) of each group is selected and compared (b).Axis in ternary plot (b) representing relative proportions of MaxRIMP for three types of variables.

Figure 7 .
Figure 7. Partial dependence plots of explanatory variables on regulating fire severity in different fire years (2000 in red and 2010 in cyan).This grouping of variables was following the rank as shown in y-axis of Figure 5.The curves demonstrate the smoothed means (solid line) and 95% confidence intervals (gray zone) of an ensemble of 50 models.Variable abbreviations were described in Table2.

Figure 7 .
Figure 7. Partial dependence plots of explanatory variables on regulating fire severity in different fire years (2000 in red and 2010 in cyan).This grouping of variables was following the rank as shown in y-axis of Figure 5.The curves demonstrate the smoothed means (solid line) and 95% confidence intervals (gray zone) of an ensemble of 50 models.Variable abbreviations were described in Table2.

Figure 8 .
Figure 8. Photographs of two unburned stands dominated by the deciduous coniferous tree Larix gmelinii (a) and the evergreen coniferous shrub Pinus pumila (b), and two 1-year post-fire stands previously dominated by Larix gmelinii (c, surface fire) and Pinus pumila (d, canopy fire) in Huzhong Natural Reserve, China.

Figure 8 .
Figure 8. Photographs of two unburned stands dominated by the deciduous coniferous tree Larix gmelinii (a) and the evergreen coniferous shrub Pinus pumila (b), and two 1-year post-fire stands previously dominated by Larix gmelinii (c, surface fire) and Pinus pumila (d, canopy fire) in Huzhong Natural Reserve, China.

Figure 9 .
Figure 9.An overview of the pre-fire vegetation composition within fires of 2000 (a) and 2010 (b).Variable abbreviations were described in Table2.

Figure 9 .
Figure 9.An overview of the pre-fire vegetation composition within fires of 2000 (a) and 2010 (b).Variable abbreviations were described in Table2.

Figure A2 .
Figure A2.The scatter plot representing inconsistent surface moisture availability (SMA) values of the whole study area (black dot) and burned pixels of 2010 fires (red dot) derived from MOD16A2 Verison 5 (V5) and Version 6 (V6) with good quality.

Figure A3 .
Figure A3.Comparison of pixel proportions of low, moderate and high severity levels between observed fire severity maps and predicted severity maps.

Figure A4 .
Figure A4.The relative importance (RIMP) of variables generated by 50 boosting regression tree (BRT) models using two different versions (V5 and V6) of surface moisture availability (SMA) derived from MODIS ET products.Variable abbreviations were described in Table2.

Figure A2 .
Figure A2.The scatter plot representing inconsistent surface moisture availability (SMA) values of the whole study area (black dot) and burned pixels of 2010 fires (red dot) derived from MOD16A2 Verison 5 (V5) and Version 6 (V6) with good quality.

Forests 2018, 9 , 26 Figure A2 .
Figure A2.The scatter plot representing inconsistent surface moisture availability (SMA) values of the whole study area (black dot) and burned pixels of 2010 fires (red dot) derived from MOD16A2 Verison 5 (V5) and Version 6 (V6) with good quality.

Figure A3 .
Figure A3.Comparison of pixel proportions of low, moderate and high severity levels between observed fire severity maps and predicted severity maps.

Figure A4 .
Figure A4.The relative importance (RIMP) of variables generated by 50 boosting regression tree (BRT) models using two different versions (V5 and V6) of surface moisture availability (SMA) derived from MODIS ET products.Variable abbreviations were described in Table2.

Figure A3 .
Figure A3.Comparison of pixel proportions of low, moderate and high severity levels between observed fire severity maps and predicted severity maps.

Forests 2018, 9 , 26 Figure A2 .
Figure A2.The scatter plot representing inconsistent surface moisture availability (SMA) values of the whole study area (black dot) and burned pixels of 2010 fires (red dot) derived from MOD16A2 Verison 5 (V5) and Version 6 (V6) with good quality.

Figure A3 .
Figure A3.Comparison of pixel proportions of low, moderate and high severity levels between observed fire severity maps and predicted severity maps.

Figure A4 .
Figure A4.The relative importance (RIMP) of variables generated by 50 boosting regression tree (BRT) models using two different versions (V5 and V6) of surface moisture availability (SMA) derived from MODIS ET products.Variable abbreviations were described in Table2.

Figure A4 .
Figure A4.The relative importance (RIMP) of variables generated by 50 boosting regression tree (BRT) models using two different versions (V5 and V6) of surface moisture availability (SMA) derived from MODIS ET products.Variable abbreviations were described in Table2.

Table 1 .
Detailed information of 21 fires included in this study.
Occurrence Date DOY † Duration (Day) Longitude Latitude Burned Area (ha)

Table 1 .
Detailed information of 21 fires included in this study.

Table 2 .
Abbreviations and descriptions of explanatory variables in this study.

Table 3 .
Accuracy assessment of predicted fire severity classification for fire year 2000 and 2010, respectively.

Table 3 .
Accuracy assessment of predicted fire severity classification for fire year 2000 and 2010, respectively.