Complex Relationships of the Effects of Topographic Characteristics and Susceptible Tree Cover on Burn Severity

Forest fires and burn severity mosaics have profound impacts on the post-fire dynamics and complexity of forest ecosystems. Numerous studies have investigated the relationship between topographic variables and susceptible tree covers with regard to burn severity. However, these relationships have not been fully elucidated, because most studies have assumed linearity in these relationships. Therefore, we examined the linearity and the nonlinearity in the relationships between topographic variables and susceptible tree covers with burn severity by comparing linear and nonlinear models. The site of the Samcheok fire, the largest recorded forest fire in Korea, was used as the study area. We generated 802 grid cells with a 500-m resolution that encompassed the entire study area and collected a dataset that included the topographic variables and percentage of red pine trees, which are the most susceptible tree cover types in Korea. We used conventional linear models and generalized additive models to estimate the linear and the nonlinear models based on topographic variables and Japanese red pine trees. The results revealed that the percentage of red pine trees had linear effects on burn severity, reinforcing the importance of silviculture and forest management to lower burn severity. Meanwhile, the topographic variables had nonlinear effects on burn severity. Among the topographic variables, elevation had the strongest nonlinear effect on burn severity, possibly by overriding the effects of susceptible fuels over elevation effects or due to the nonlinear effects of topographic characteristics on pre-fire fuel conditions, including the spatial distribution and availability of susceptible tree cover. To validate and generalize the nonlinear effects of elevation and other topographic variables, additional research is required at different fire sites with different tree cover types in different geographic locations.


Introduction
Forest fires significantly impact forest ecosystems by creating unique burn mosaics and are considered key components of the spatiotemporal dynamics and complexity of forest ecosystems [1][2][3][4]. Forest fires are significantly associated with landslides [5], soil microbial processes [6], seed germination and sprouting [7], live and dead vegetation structure [8,9] and composition [10] and carbon budgets [11]. In particular, burn severity is a critical factor in understanding the degree of influence of forest fires on forest ecosystems, post-fire vegetation responses and heterogeneity of vegetation composition and configuration in burned areas by affecting the availability of seed sources, sprouting rates, soil humidity, soil nutrients, lights, wind speed, alien plant invasion, tree mortality, animal populations and community dynamics [10,[12][13][14][15] at various spatial and temporal scales. Therefore, burn severity is the most critical factor in determining the dynamics and complexity of the post-fire response of a damaged forest ecosystem to a fire event [10,16]. Due to the strong impacts of burn severity on forest ecosystems, understanding how burn severity is influenced by various environmental factors during a fire can offer profound insights to minimize economic loss, preserve forest ecosystems, practice fire-resilient forest management strategies and plan effective post-fire restoration for planners and decision makers. Although numerous studies have been conducted, the relationships between burn severity and environmental factors remain unknown due to the dynamics and complexity of these relationships and interactions among environmental factors during fire events. Considering that recent trends in warming and increasing aridity could increase the probability of weather associated with severe fires [17], there is increasing concern over the management of burn severity and environmental factors escalating fire risks.
The three primary factors (i.e., the fire triangle) affecting burn severity are topography (e.g., elevation, slope, aspect, topographic position, solar radiation and topographic wetness), weather (e.g., temperature, relative moisture, wind speed, drought factor and regional climate) and fuel (e.g., forest cover type, composition and configuration of forest, shrubs, heterogeneity and tree density) [10,18]. Some studies have shown that topographic characteristics are the most significant variables affecting burn severity. For example, Fang et al. [19] reported that topography was one of the most significant factors explaining burn severity in the boreal forest of the Great Xing'an Mountains, China. Similarly, Bigler et al. [20] found elevation to be the most significant predictor of fire severity in Rocky Mountain subalpine forests. To understand the effects of weather on burn severity, two different types of weather should be considered: climate and fire weather. Climate reflects broad weather patterns and is closely associated with long-term annual and seasonal variation in precipitation, temperature, relative moisture, solar radiation, drought and wind on a regional scale [19,21]. In general, higher temperatures, lower humidity and dry conditions create environments that are more susceptible to forest fires due to decreased fuel moisture and the lower energy required for pre-heating [22]. Meanwhile, fire weather reflects site-specific weather conditions during burning incidents and there is an abundance of evidence indicating that fire weather conditions (i.e., temperature, relative moisture, wind speed and drought factor) are the most critical factors in explaining many fire characteristics such as ignition, spread and size and burn severity [20,23,24]. Based on fire studies, forest fires are highly selective [3]. For example, less susceptible forest cover types (e.g., olive groves) can constrain fire spread and lower the burn severity, whereas susceptible forest covers (e.g., aspen, red pine and birch) can enhance the spread of fire and burn severity [3]. Moreover, many fuel traits directly impact burn severity. Some of the primary traits of fuels related to burn severity are spatial distribution [10], density [25][26][27], moisture [28], fuel type [3,29] and pre-fire forest heterogeneity [10].
Reviewing the literature on fire with regard to variables affecting burn severity reveals the complexity of these relationships. Despite the large body of research, the influence of the fire triangle on burn severity remains poorly understood due to the complex nature of the relationships. In the literature, some of this complexity is due to interactions among topography, weather and fuel and their spatially varying effects [30]. Simultaneously, this complexity is partially due to potential nonlinear relationships between burn severity with topography, weather and fuels [3,19].
From this perspective, we examined the presence of nonlinear relationships of topographic characteristics and fuels with burn severity. Several studies have suggested the possibility of nonlinear relationships between burn severity and topography, weather and fuels. For example, Lee et al. [3] reported the possibility of nonlinear relationships between burn severity and elevation and slope in the hilly forest of Samcheok, Korea. Recently, Fang et al. [19] also reported nonlinear response curves of the top three explanatory variables for both fire size and fire severity models in the boreal forest of the Great Xing'an Mountains, China. The underlying rationale of this study was that the complexity and dynamic effects of the fire triangle were in part due to the presence of nonlinear relationships. In this study, we focused on topographic characteristics and susceptible tree cover type. Weather conditions were not considered because we examined only one study site and all analysis units shared the same climate conditions; additionally, the local fire weather conditions within the study areas were not available. Despite these limitations, the results of this study expand our knowledge of the relationships of burn severity with topography and fuels and help to clarify the nature of these relationships, which is crucial for informing forest managers and policy makers.

Study Site Characteristics and Samcheok Forest Fire
Samcheok is located in the eastern region of the Korean Peninsula and borders the East Sea. The study area experiences a typical monsoon climate, characterized by substantial variation in precipitation and temperature. In particular, precipitation is concentrated (about two-thirds of the annual total) in summer and there is a dry season in spring and fall. Overall, the average annual temperature and precipitation are 12.1 • C and 1342 mm, respectively [30]. The fire history indicates that 75% of forest fires in Korea occur in spring (March-May) due to dry and windy weather conditions [31]. Several studies have documented the complex topography of the site (e.g., [3,30]). Before the fire, 65% of the site was dominated by Japanese red pine (Pinus densiflora), which is the most susceptible forest cover type in Korea. The forest had an average age of approximately 30 years and was at the secondary succession stage, following the Korean War in the early 1950s [3]. The Samcheok fire occurred on 7 April 2000, starting from a garbage burn site and burned for 9 days. It is the largest forest fire recorded in Korea since 1960 [32]. During the fire event, the maximum wind speed was 26.8 m/s and the minimum relative humidity was 7%, representing fire-prone conditions. Such susceptible forest cover type and fire-prone weather conditions resulted in 16,151 ha of unique burned mosaic in one of the most preserved landscapes in Korea ( Figure 1). Several government agencies including the Korea Fire Service (KFS) [31], National Institute of Forest Science, Ministry of Environment and a number of private institutions have set long-term monitoring sites in the area to monitor post-fire changes in various biotic and abiotic factors such as vegetation, soil, stream water quality, insects, butterflies, birds and mammals.

Spatial Unit of Analysis and Mapping the Percentage of Red Pine Trees
To capture the spatial distribution of Japanese red pine trees, topographic characteristics (i.e., elevation, slope, topographic wetness index [TWI] and solar radiation index [SRI]) and burn severity in the study area, we divided the study area into 500-m grid cells. We generated 802 grid cells (500 m) covering the entire study areas in the computer-aided draft (CAD) of Autodesk and converted the grid CAD file into a polygon shape file using the geographic information system (GIS). A unique identification code was assigned to each grid cell and the mean value of all selected variables within each grid cell was computed with GIS using the overlay function. Previous studies have shown that this spatial resolution can effectively capture complex topographic characteristics and burn severity (e.g., [10,30]). Forest fire is very selective to forest cover type. Japanese red pine is the most susceptible forest cover type in Korea and other forest cover types including mixed forest or broad-leaved forest have very weak relationships with burn severity [3,30]. The high susceptibility of red pine trees has been reported in other geographic locations (e.g., [20,33]). Red pine is the dominant forest cover type in Korea. Therefore, the primary fire suppression strategies of the KFS have involved the spatiotemporal management of red pine and we only considered red pine as the fuel type in explaining the variance of burn severity. We used the National Forest Classification Digital Map (NFCDM) released in 2000 by the Korea Forest Service (KFS) to compute the percentage of red pine trees within each grid cell. The NFCDM was originally a polygon shape file, which we converted into

Spatial Unit of Analysis and Mapping the Percentage of Red Pine Trees
To capture the spatial distribution of Japanese red pine trees, topographic characteristics (i.e., elevation, slope, topographic wetness index [TWI] and solar radiation index [SRI]) and burn severity in the study area, we divided the study area into 500-m grid cells. We generated 802 grid cells (500 m) covering the entire study areas in the computer-aided draft (CAD) of Autodesk and converted the grid CAD file into a polygon shape file using the geographic information system (GIS). A unique identification code was assigned to each grid cell and the mean value of all selected variables within each grid cell was computed with GIS using the overlay function. Previous studies have shown that this spatial resolution can effectively capture complex topographic characteristics and burn severity (e.g., [10,30]). Forest fire is very selective to forest cover type. Japanese red pine is the most susceptible forest cover type in Korea and other forest cover types including mixed forest or broad-leaved forest have very weak relationships with burn severity [3,30]. The high susceptibility of red pine trees has been reported in other geographic locations (e.g., [20,33]). Red pine is the dominant forest cover type in Korea. Therefore, the primary fire suppression strategies of the KFS have involved the spatiotemporal management of red pine and we only considered red pine as the fuel type in explaining the variance of burn severity. We used the National Forest Classification Digital Map (NFCDM) released in 2000 by the Korea Forest Service (KFS) to compute the percentage of red pine trees within each grid cell. The NFCDM was originally a polygon shape file, which we converted into GRID file format with a 50-m resolution, from which the percentage of red pine trees within each grid cell was computed.

Mapping Topographic Characteristics
The topographic characteristics of the study site are very complex at small and intermediate scales ( Figure 1). In fire research, various indicators have been adopted to delineate topographic characteristics. The most common topographic indicators used in studies include elevation, slope, aspect, topographic wetness index (TWI), solar radiation index (SRI), topographic position index, elevation relief ratio, heat load index, topographic roughness index and gullies (e.g., [19,[34][35][36][37][38][39][40][41]). Because there is no clear evidence indicating which variables are more effective in delineating topographic characteristics, we tried to adopt indicators that have proven to be effective in capturing the complex topographic characteristics of our study areas in previous studies (e.g., [3,10,30]). From an implementation perspective, simple indicators are better for forest managers and decision makers who might be unfamiliar with complex topographic indicators. At the same time, most topographic indicators are strongly correlated to each other, resulting in multicollinearity and using many indicators does not necessarily enhance estimated model performance. Therefore, we focused on major topographic characteristics including altitude, steepness, topographic type (i.e., valley/ridge type) and the effects of incident solar radiation on the surface or fuels.
Based on these considerations, the selected indicators for the study included elevation (m), slope ( • ), TWI and SRI. Elevation (m) and slope ( • ) are direct measurements that delineate topographic characteristics and have shown strong relationships with burn severity in many fire studies (e.g., [10,19,20,30,42,43]). TWI and SRI are relative measures of topographic characteristics. TWI is a relative measure of the long-term soil moisture availability based on local slope and upslope areas draining at a certain point (see [44,45]). Slope aspect is an important factor for understanding potential direct incident radiation and temperature at a given location. However, using a simple aspect is somewhat problematic because it cannot distinguish between 1 • and 360 • [38]. Therefore, the SRI was proposed to compute the potential annual direct incident solar radiation at a given point without the limitations of simple aspect (see [46]). A study on formation of the boundary of burned areas found TWI to be a significant factor in higher-elevation areas and SRI was an important factor in areas with low susceptible forest covers [10]. Therefore, TWI and SRI appear to have important roles in fire behavior under somewhat different local conditions. Most forest fire studies have used topographic analysis and computing topographic index based on digital elevation models (e.g., [19,47]). Similarly, we used a digital elevation model to analyze the mean elevation, slope, TWI and SRI within the 500-m grid cells.

Mapping Burn Severity
Even though the term fire severity has been used widely and interchangeably with burn severity in the literature [26,48] and the two phenomena share some commonalities in fire and post-fire effects, there are critical differences between them. For example, fire severity is associated with active fire characteristics during the fire event and the immediate post-fire effects on the local environment. Meanwhile, burn severity incorporates both short-and long-term post-fire effects on local and regional environments [27]. In this study, burn severity referred to the degree of environmental change (i.e., change in vegetation reflectance) caused by fire [26,30,48], allowing burn severity to be measured based on a number of variables. One of the most common variables used to measure burn severity is the vegetation consumed by fire in the damaged area, which we adopted for burn severity in this study. Burn severity refers to the degree of change in vegetation reflectance captured in pre-and post-fire satellite images in this study. The difference in the near-infrared reflectance of satellite images taken pre-and post-fire is the delta normalized burn ratio (dNBR), which has been used to capture burn severity effectively in many forest fire studies (e.g., [9,33,49]). A number of studies have shown that dNBR very effectively captures burn severity in the same study areas and have reported a strong association between dNBR and the field-measured composite burn index (e.g., [10,30]). We used Landsat Thematic Mapper satellite images with a 30-m resolution taken on 7 May 1999 (pre-fire) and 5 May 2000 (post fire) to compute the dNBR and the satellite-derived burn severity was classified into six classes (extreme severity, very high severity, high severity, moderate severity, low severity and unburned) using the Remote Sensing Application Center [49] criteria. Then, burn severity derived from satellite images was validated with field-measured semi-composite burn index values (r = 0.832, p < 0.01 [3]). Finally, we computed the mean burn severity within each 500-m grid cell for the analysis.

Linear and Nonlinear Model Estimation
In forest fire studies, simple correlations and linear regression models are the most common conventional approaches to investigate the relationships between burn severity and environmental variables, relying on the underlying assumption of a linear relationship. However, this assumption is highly questionable [10,50]. Linear correlation and regression are useful when quantifying the magnitude, direction and significance of the relationships between burn severity and environmental variables such as topographic variables and fuel types. However, conventional linear approaches may not accurately represent the true nature of the relationships between burn severity and environmental variables, which could mislead forest managers and policy makers. To avoid nonlinearity issues in dealing with the fire regime and environmental variables, a few approaches have been proposed such as artificial neural networks [51], stochastic processing [52] and spatial clustering [53]. Recently, flexible regression models, such as generalized additive models (GAMs; [54]), have been proposed to handle the nonlinear effects of continuous covariates on the response variable [50,55].
We used the R-package (R Core Development Team) to estimate the linear and nonlinear models of the relationship between burn severity with topographic variables and red pine trees. We adopted the GAM in the R-package to fit nonlinear models, as GAM is very flexible and can provide an excellent fit for nonlinear relationships [56]. It is a multiple regression model in which the additive nature of the model is maintained but the linear regression lines are replaced by nonparametric functional curves with multiple parameters. An interesting feature of GAM is that it allows determination of the shape of the response curves from the data instead of fitting an a priori parametric model, making it data-driven instead of model driven [57]. GAMs have been widely applied to investigate nonlinear relationships between dependent and non-explanatory variables in various fields of study such as plant and aquatic ecology (e.g., [58]) and water quality (e.g., [59]), as well as for nonlinear relationships among biotic and bon-biotic variables in environmental settings due to complex interactions between environmental factors. Many studies have identified nonlinear relationships among components of forest ecosystems (e.g., [3,50,[60][61][62][63][64][65]). In particular, Lee et al. [3] reported that burn severity might be related to topographic characteristics such as elevation and slope. In their study, nonlinear models performed much better than linear models in explaining the variances of burn severity with elevation and slope due to the spatial variances in fuel moisture, temperature, precipitation and canopy bulk density [66]; the compounding effects of fuels and the chimney effect [22]; and water potential [67].
The general linear multiple regression model (LM) for burn severity with topographic variables and susceptible fuel is: where LM bs is burn severity; α is the intercept; β i is the coefficient of elevation, slope, TWI, SRI and the percentage of red pine trees; and ε is the error term. In this equation, β i represents the linear contributions of elevation, slope, TWI, SRI and the percentage of red pine trees to burn severity. However, the GAM assumes that the contribution of each variable varies over the range of each variable, which can be represented by a function instead a single coefficient. Applying this assumption to Equation (1) yields the following conceptual equation of GAM for burn severity: where g is the unknown identity link; GAM bs is the burn severity; α is the overall intercept; f i is the unknown smooth functions of elevation, slope, TWI, SRI and the percentage of red pine trees; and ε is the error term. In this context, the GAM for burn severity is the sum of the functions of elevation, slope, TWI, SRI and the percentage of red pine trees within each grid cell. Therefore, each variable can have varying coefficients defined as a function and the contributions of each variable to burn severity can differ among grid cells depending on the settings of the topographic characteristics (i.e., elevation, slope, TWI and SRI) and susceptible fuel (i.e., percentage of red pine trees). We estimated the LM and nonlinear GAM using the R-package to explain the variances in burn severity with topographic characteristics (i.e., elevation, slope, TWI and SRI) and susceptible fuels type (i.e., percentage of Japanese red pine trees). Then we compared the performance of the two estimated models using the coefficient of determination (R 2 ), Akaike's information criterion (AIC) and the Bayesian information criterion (BIC). In general, models with higher R 2 values explain the variance of the explanatory variables better than models with lower R 2 values. Similarly, the AIC represents the relative likelihood of delineating the true nature of the relationships among variables and smaller AIC values represent a higher likelihood of explaining the nature of the data than other models with higher AIC values; therefore, larger AIC values for models are associated with a lower probability [59]. Despite its usefulness and widespread use, AIC has been criticized for its tendency to select overly complex models [68]. In addition, AIC is inconsistent and neglects the sampling variability of the parameters [69]. As an alternative, BIC can be used as an indicator of model selection without these problems [68][69][70]. The BIC value of a model is an estimate of a function of the posterior probability that the model is true; thus, lower BIC values indicate that a model is more likely to be the true model. Despite their different theoretical backgrounds, benefits and drawbacks, lower AIC and BIC values both indicate better model performance in a practical sense.

Descriptive Statistics and Spatial Distribution
The distribution of burn severity at the Samcheok site followed a near-normal distribution, with a large number of observations around the mean value (mean = 3.88, standard deviation = 0.66). Only 2 grid cells experienced low burn severity (<2), whereas 26 grid cells experienced very high severity and extremely severity. The majority of grid cells were categorized as high severity, moderate severity and low severity. Previous studies have documented the complex topography of the study area (e.g., [3,10,30]). The maximum and minimum elevations of the study areas were 849.5 m and 0.0 m, respectively. Some parts of the damaged area shared a boundary with the East Sea and there was great variance in elevation (Figures 1 and 2). However, the distribution of elevation and slope differed substantially. The observed elevation showed a left-skewed distribution (skewness = 0.825, kurtosis = −0.25), whereas the slope showed a somewhat normal distribution (skewness = −0.235, kurtosis = −0.421). The TWI and SRI showed somewhat opposing distribution shapes, in which TWI had a left-skewed distribution (skewness = 0.833, kurtosis = 1.25) and SRI had a right-skewed distribution (skewness = −0.53, kurtosis = 0.03). The mean TWI and SRI values were 2.45 and 0.88, respectively. Interestingly, the percentage of red pine trees within the grid cells showed a very different distribution than the other variables ( Figure 3). In total, 61 grid cells (7.8%) had no pine trees, whereas 77 grid cells were entirely covered with red pine trees (9.9%). The mean percentage of red pine trees within a grid was 59.02, which suggests a high pre-fire dominance of red pine trees in the study area. . Burn severity, slope, TWI and SRI exhibited near-normal distributions, whereas elevation was left-skewed. Interestingly, the distribution of the percentage of red pine trees showed that the frequency of extreme values at the minimum (0%) and maximum (100%) values were very high, indicating high variance of red pine trees over the study area.

Correlation Analysis and Scatter Plots
The results of the correlation analyses and scatter plots indicated a strong association of burn severity with topographic characteristics and percentage of red pine trees. Burn severity was negatively correlated with elevation (r = −0.25, p < 0.01) and TWI (r = −0.15, p < 0.01) and positively correlated with the percentage of red pine trees (r = 0.47, p < 0.01) ( Table 1). Burn severity was not significantly correlated with slope or SRI. However, the scatter plots in Figure 3 were strongly suggestive of nonlinear relationships between burn severity and elevation and slope. The scatter plot between burn severity and the percentage of red pine trees seemed to show a linear relationship. Correlation analysis indicated that burn severity was higher in low areas covered with more red pine trees. Elevation and slope showed a high positive correlation (r = 0.69, p < 0.01) but the shape of the scatter plot suggested the presence of a nonlinear relationship between elevation and slope. Elevation showed a strong negative correlation with the percentage of red pine trees (r = −0.42, p < 0.01), which suggests that red pine trees were more common in lowland areas. However, it was unclear whether the relationship was linear or nonlinear based on the shape of the scatter plot ( Figure 4). Simultaneously, the percentage of red pine trees showed a negative relationship with slope (r = −0.12, p < 0.01), which suggests that mild slopes were likely covered with more red pine trees than steep slopes. There was a strong negative correlation between elevation and TWI (r = −0.57, p < 0.01), which suggests that lowlands likely had higher wetness values. The scatterplot between elevation and TWI appeared as a gently curved shape, which suggests a nonlinear relationship between elevation and TWI. In addition, TWI showed a strong negative correlation with slope (r = −0.65, p < 0.01), which indicates that areas with steep slopes might have low TWI values. The shape of the scatter plot between TWI and slope appeared to be relatively linear.

Correlation Analysis and Scatter Plots
The results of the correlation analyses and scatter plots indicated a strong association of burn severity with topographic characteristics and percentage of red pine trees. Burn severity was negatively correlated with elevation (r = −0.25, p < 0.01) and TWI (r = −0.15, p < 0.01) and positively correlated with the percentage of red pine trees (r = 0.47, p < 0.01) ( Table 1). Burn severity was not significantly correlated with slope or SRI. However, the scatter plots in Figure 3 were strongly suggestive of nonlinear relationships between burn severity and elevation and slope. The scatter plot between burn severity and the percentage of red pine trees seemed to show a linear relationship. Correlation analysis indicated that burn severity was higher in low areas covered with more red pine trees. Elevation and slope showed a high positive correlation (r = 0.69, p < 0.01) but the shape of the scatter plot suggested the presence of a nonlinear relationship between elevation and slope. Elevation showed a strong negative correlation with the percentage of red pine trees (r = −0.42, p < 0.01), which suggests that red pine trees were more common in lowland areas. However, it was unclear whether the relationship was linear or nonlinear based on the shape of the scatter plot ( Figure 4). Simultaneously, the percentage of red pine trees showed a negative relationship with slope (r = −0.12, p < 0.01), which suggests that mild slopes were likely covered with more red pine trees than steep slopes. There was a strong negative correlation between elevation and TWI (r = −0.57, p < 0.01), which suggests that lowlands likely had higher wetness values. The scatterplot between elevation and TWI appeared as a gently curved shape, which suggests a nonlinear relationship between elevation and TWI. In addition, TWI showed a strong negative correlation with slope (r = −0.65, p < 0.01), which indicates that areas with steep slopes might have low TWI values. The shape of the scatter plot between TWI and slope appeared to be relatively linear. Overall, the results presented in Table 1 suggested that the topographic characteristics, percentage of red pine trees and burn severity were closely associated. The percentage of red pine trees was closely associated with elevation and TWI showed a strong association with elevation and slope, suggestive of multicollinearity among elevation, slope and TWI. The comparison of the correlations and scatter plots provided interesting insights into the relationships of burn severity with topographic characteristics and red pine trees. The scatter plots in Figure 3 suggested the high possibility of the presence of a nonlinear relationship between burn severity and some of the observed variables, which did not appear to be significant in the correlation analysis. Because correlation analyses delineate linear relationships among variables, simple correlation analyses were not sufficient to capture the true nature of burn severity with the topographic characteristics and percentage of red pine trees at the Samcheok fire site.  Overall, the results presented in Table 1 suggested that the topographic characteristics, percentage of red pine trees and burn severity were closely associated. The percentage of red pine trees was closely associated with elevation and TWI showed a strong association with elevation and slope, suggestive of multicollinearity among elevation, slope and TWI. The comparison of the correlations and scatter plots provided interesting insights into the relationships of burn severity with topographic characteristics and red pine trees. The scatter plots in Figure 3 suggested the high possibility of the presence of a nonlinear relationship between burn severity and some of the observed variables, which did not appear to be significant in the correlation analysis. Because correlation analyses delineate linear relationships among variables, simple correlation analyses were not sufficient to capture the true nature of burn severity with the topographic characteristics and percentage of red pine trees at the Samcheok fire site.

Estimating Simple LMs and GAMs
To examine the causal relationships of the topographic variables and percentage of red pine trees with burn severity, we estimated LMs and GAMs for each non-explanatory variable (i.e., elevation, slope, TWI, SRI and pine trees) separately. Simple LMs and GAMs of each variable were compared using the R 2 , AIC and BIC values to determine whether nonlinear relationships were present between each variable and burn severity. If the GAM outperformed the LM in explaining burn severity, the presence of a nonlinear relationship between non-explanatory variables (i.e., topographic variables and percentage of red pine trees) and the explanatory variable (i.e., burn severity) was assumed. Table 2 summarizes the estimated simple LMs and GAMs. The results strongly indicated that the GAMs of all topographic variables outperformed the LMs in terms of R 2 , AIC and BIC values. In particular, the GAM (R 2 = 0.23, AIC = 1,427.9, BIC = 1,467.7) of elevation performed much better than the LM (R 2 = 0.06, AIC = 1581.1, BIC = 1595.5). Similarly, the GAM (R 2 = 0.10, AIC = 1550.5, BIC = 1598.4) of slope showed a better performance than the LM (R 2 = 0.00, AIC = 1630.8, BIC = 1644.8). The GAMs of TWI and SRI outperformed the respective LMs. However, the low R 2 values of both the GAMs and LMs of slope, TWI and SRI indicated that their influences on burn severity were modest at best. Among the topographic variables, elevation appeared to have the highest influence on burn severity. However, the influence of elevation on burn severity was not constant and varied over the elevation range (i.e., it was nonlinear). Interestingly, the R 2 , AIC and BIC values of the LM and GAM of the percentage of red pine trees were almost identical, which did not support the presence of a nonlinear relationship. In summary, the percentage of red pine trees was linearly associated with burn severity, while elevation showed a significant nonlinear relationship with burn severity. Slope, TWI and SRI also showed modest nonlinear relationships with burn severity. Table 2. Estimated simple LMs and GAMs. The low F-values of the estimated LMs of slope and SRI indicated that the estimated simple LM models were not significant. No differences in the effects of red pine trees on burn severity were observed between the LM and GAM. The estimated models of red pine trees were almost identical in terms of F-value, adjusted R 2 , AIC and BIC.  Figure 4 presents a graphical representation of the smooth functions of each topographic variable. Elevation and TWI had relatively simple nonlinear relationships with burn severity, whereas slope and SRI have relatively complex nonlinear relationships with burn severity. The smooth function of elevation (Figure 4) could be classified into two rough regions: a positive influence region and negative influence region. In Region 1, low elevations showed escalated burn severity and the steep positive function line (i.e., solid line in the middle of the shaded area) suggested that elevation had strong positive impacts on burn severity within the region. However, the influence of elevation on burn severity in Region 2 gradually changed towards a negative influence, which suggests that burn severity decreased with increasing elevation. Compared to elevation, the slope showed a more complex smooth function, which changed directions multiple times (Figure 4). Slope lowered burn severity in Regions 1 and 3 but increased burn severity in Regions 2 and 4. However, only a small number of observations were found in Regions 1 and 4 and their confidence intervals were relatively wide, suggesting that the direction of the influence of slope on burn severity might not be highly significant.
The smooth function of TWI appeared to show a relatively simple, positive influence at low TWI values (Region 1) and a negative influence at high TWI values (Region 2). The smooth function of SRI could be classified into three regions: Region 1 (negative influence), Region 2 (positive influence) and Region 3 (negative influence). Notably, Region 1 (negative effect) and Region 4 (positive effect) in the smooth function of slope and Region 1 (negative effect) in the smooth function of SRI were unlikely to be overestimated, because there were not many cases and their R 2 values were extremely low (Table 3). Overall, the smooth functions of the topographic variables revealed inconsistent influences of topographic characteristics on burn severity, which varied by elevation, slope, TWI and SRI. Elevation showed the strongest nonlinear effects on burn severity. However, the effects of slope, TWI and SRI on burn severity were modest at best. By contrast, the susceptible forest cover type (i.e., red pine trees) showed a linear relationship with burn severity. This emphasizes the importance of controlling the density and amount of flammable forest cover to minimize burn severity during fire events and manage fire-resilient forests. of SRI were unlikely to be overestimated, because there were not many cases and their R 2 values were extremely low (Table 3). Overall, the smooth functions of the topographic variables revealed inconsistent influences of topographic characteristics on burn severity, which varied by elevation, slope, TWI and SRI. Elevation showed the strongest nonlinear effects on burn severity. However, the effects of slope, TWI and SRI on burn severity were modest at best. By contrast, the susceptible forest cover type (i.e., red pine trees) showed a linear relationship with burn severity. This emphasizes the importance of controlling the density and amount of flammable forest cover to minimize burn severity during fire events and manage fire-resilient forests.

Estimating Multiple LM and GAM
We estimated multiple LM and GAM with all of the independent variables (i.e., percentage of red pine trees, elevation, slope and TWI) to explain the variance of burn severity using the R-package. Multiple LM based on the least square method can provide the mean effects of an independent variable on explanatory variables while controlling for other independent variables in the model. Unlike LMs, GAMs do not provide coefficients for each independent variable because a single coefficient is pointless if an independent variable is related to explanatory variables in a nonlinear fashion. Instead, they provide an intercept and the approximate significance of smooth terms (i.e., significance of the nonlinear relationship, (see [54] for details) of each independent variable with explanatory variables and deviance explained by the model.
The best fitting LM was selected using a stepwise method based on the R 2 , AIC and BIC values. The estimated best-fit model included elevation, TWI and the percentage of red pine trees, while slope and SRI did not have significant effects on burn severity ( Table 3). The estimated multiple linear model indicated that the percentage of red pine trees (b = 0.007, β = 0.34, p < 0.01) escalated burn

Estimating Multiple LM and GAM
We estimated multiple LM and GAM with all of the independent variables (i.e., percentage of red pine trees, elevation, slope and TWI) to explain the variance of burn severity using the R-package. Multiple LM based on the least square method can provide the mean effects of an independent variable on explanatory variables while controlling for other independent variables in the model. Unlike LMs, GAMs do not provide coefficients for each independent variable because a single coefficient is pointless if an independent variable is related to explanatory variables in a nonlinear fashion. Instead, they provide an intercept and the approximate significance of smooth terms (i.e., significance of the nonlinear relationship, (see [54] for details) of each independent variable with explanatory variables and deviance explained by the model. The best fitting LM was selected using a stepwise method based on the R 2 , AIC and BIC values. The estimated best-fit model included elevation, TWI and the percentage of red pine trees, while slope and SRI did not have significant effects on burn severity ( Table 3). The estimated multiple linear model indicated that the percentage of red pine trees (b = 0.007, β = 0.34, p < 0.01) escalated burn severity while elevation (b = −0.001, β = −0.27, p < 0.01) and TWI (b = −0.305, β = −0.29, p < 0.01) decreased burn severity. The F-value of the estimated the best-fit model was 97.63 (p < 0.01) and the low variance inflation factor (VIF) values of all of the independent variables indicated a lack of multicollinearity among the variables included in the model. The estimated multiple model explained about 27% of the variance of burn severity and a large proportion of the variance of burn severity was not explained by these three independent variables. The estimated multiple best-fit LM was: LM bs = 4.47 + (−0.001 × elevation) + (−0.31 × TWI) + (0.007 × percentage of red pine) In estimated GAMs, the nonlinearity between an independent variable and explanatory variables can be judged with the expected default frequency (EDF) value of each variable. When the EDF equals 1, the explanatory variables and independent variable have linear relationships, whereas higher EDF values are an indication of increasingly nonlinear relationships [54]. The estimated GAM in Table 4 (family = Gaussian, link function = identity) indicated that all of the topographic characteristics were related with burn severity in a nonlinear form, whereas the percentage of red pine trees (EDF = 1) was associated with burn severity in a linear form. TWI (8.33, p < 0.001) and elevation (6.23, p < 0.001) showed relatively higher EDF values, suggestive of a strong nonlinear relationship with burn severity. Among the topographic variables, the slope appeared to have the lowest EDF value (1.95, p < 0.001), which indicates a relatively less nonlinear relationship with burn severity. The EDF of the percentage of red pine trees was 1, which indicates a linear relationship with burn severity. The estimated GAM for burn severity was: g (GAM bs ) = 3.64 + (0.004 × percentage of red pine) + s (elevation) + s (slope) + s (TWI) + s (SRI) (4) where g is the link function (identity) and s is the smooth function of each variable. The adjusted R 2 of the estimated GAM for burn severity was 0.36 (Table 4), which strongly suggested that the GAM outperformed the LM (adjusted R 2 = 0.27, Table 2) in explaining the variance of burn severity at the Samcheok burn site. In addition, the AIC value of the GAM was 1297.43, which was considerably smaller than that of the LM (1385.64). The BIC value (1402.22) of the GAM was also smaller than that of the LM (BIC = 1422.23), indicating that the GAM performed better than the LM in explaining burn severity. Analysis of variance (ANOVA) was performed to test for significant differences between the two models in explaining burn severity. The results indicated that there was a significant difference between the two models in explaining burn severity ( Table 5). The ANOVA and comparisons of the R 2 , AIC and BIC of the two models strongly indicated that the GAM performed better than the LM in explaining burn severity; therefore, topographic characteristics such as elevation, slope, TWI and SRI had nonlinear relationships with burn severity. However, the percentage of red pine trees appeared to have a linear relationship with burn severity, which implies that more pine trees resulted in greater burn severity regardless of the topographic characteristics. Table 4. Estimated GAM. About 36% of burn severity could be explained by the estimated GAM. The percentage of red pine trees showed a linear relationship with burn severity, whereas all of the other independent variables exhibited non-linear relationships with burn severity.

Linear Relationship of Red Pine Trees with Burn Severity
The results of this study strongly support the findings of previous studies reporting close associations between red pine trees and burn severity. All of the analysis results in this study including the correlation analysis (Table 1), estimated LM (Table 2) and GAM ( Table 3) strongly suggested that higher percentages of red pine trees intensified the degree of burn severity, in agreement with numerous previous studies (e.g., [3,10,18,26]).
In particular, a comparison of the β-values of independent variables in the estimated LM highlighted that the percentage of red pine trees was the primary determinant of burn severity. A few studies also reported a stronger association between burn severity and pre-fire forest cover type than topography or fire weather (e.g., [19,20]). For example, Lee et al. [3] reported that the percentage of red pine trees was the most significant factor in explaining burn severity compared to topographic characteristics (i.e., elevation and slope) using regression tree analysis and landscape metrics. In their regression tree, the percentage of red pine trees appeared as the first-order variable in explaining burn severity, whereas the spatial configuration of red pine trees and topographic variables were the second and third order variables, respectively. Interestingly, Fang et al. [19] quantified the relative importance of fire weather, fuel type and topography at the Great Xing'an Mountain fire site in northeast China. The authors reported that the spatial extent of burned area was determined primarily by the weather conditions during the fire event. However, burn severity was mostly influenced by fuel type and the weather condition during the fire event was a relatively less significant factor of burn severity than forest cover type and topography. Bigler et al. [20] also reported that the most significant predictors of burn severity were pre-fire vegetation conditions and elevation in Rocky Mountain subalpine forests.
In addition, we found that the percentage of red pine trees in the grid cells was linearly associated with burn severity and a higher percentage of red pine trees in grids cells resulted in higher burn severity. The study area has complex and dynamic topographic characteristics, with great variation in elevation, slope, TWI and SRI (Figure 3). According to the KFS [31] fire report, there were substantial variations in wind speed and relative humidity during the fire event. Despite the variations in topographic characteristics and changes in weather conditions, forests with a higher percentage of red pine trees burned more severely than those with fewer pine trees.
However, we were unable to find many studies on the linear and nonlinear effects of susceptible tree cover types on burn severity. Ríos-Pena et al. [50] modeled the occurrence of wildfires in Galicia, Spain, using binary structured additive regression. Although they did not directly examine the nonlinear effects of susceptible fuels on burn severity, their results suggested a linear relationship between red pine trees and burn severity. In addition, the study suggested that different fuel types had somewhat mixed linear and nonlinear effects on the occurrence of wildfire, where all tree-based fuels, including shrubs and wood residuals, had linear effects, while litter under trees had nonlinear effects. However, they did not provide clear explanations for the differences in the effects among fuel types. Moreover, elevation showed declining negative nonlinear effects on burn severity and the nonlinear effects of litter were due to the nonlinear effects of elevation. We observed similar linear effects of susceptible forest cover types. Considering that the burn process requires consumable fuels, this was unsurprising and emphasizes the importance of silviculture and forest management for fire-resilient forests, because fuel characteristics are the only manageable factors in the forest fire triangle (i.e., fuel, topography and fire weather).
The linear relationship between susceptible forest cover and burn severity in lowering the density of flammable trees in forests should be a priority in forest management, because the degree of burn severity is expected to rapidly decrease in response to a decrease in flammable trees in forests. However, severe fire weather conditions can overwhelm fuel and topographic characteristics due to a stronger association between weather and fire mechanisms and its higher variability compared to fuels [19,71]. Therefore, a strong linear relationship between red pine trees and burn severity might only occur under moderate weather conditions.

Non-Linear Relationship of Topographic Characteristics with Burn Severity
Topographic characteristics have direct impacts on burn severity and fire behavior during fire events. However, the relationships between burn severity and topographic characteristics (e.g., elevation, slope and aspect) remain somewhat controversial. For example, Lee et al. [3] found negative correlations between burn severity and elevation in Samcheok, Korea. Similarly, Weatherspoon and Skinner [72] reported that higher elevation was associated with lower burn severity because of cooler temperatures and higher humidity. In contrast, other studies have reported a positive association of burn severity with elevation and slope [34,73,74]. For example, Fang et al. [19] reported a strong association of severe burning with high elevations and steep slopes in the boreal forest of the Great Xing'an Mountains, China. The contradictory results of the effects of topography on burn severity are in part due to the complex interactions of topography with fuels and fire weather [26,29,74], spatially varying effects [30] and indirect effects. In addition, topography and burn severity may be associated in a nonlinear manner [10]. Bigler et al. [20] reported that the local effects of topography declined with increasingly severe fire weather, particularly across short elevation gradients. Therefore, it is difficult to generalize or simplify the relationships of topography with burn severity due to contradictory evidence from studies of forest fires and the complex nature of these relationships.
The effects of topographic variables including elevation, slope, TWI and SRI on burn severity in this study showed nonlinear relationships with burn severity, even changing from negative to positive relationships in some cases. The nonlinear effects of topographic characteristics in this study indicate that their influences on burn severity are not constant but rather, significantly vary in terms of the direction and degree of their influences. For example, the influences on burn severity at one location might differ from those at different locations at the same fire site. Therefore, the influences of topographic characteristics on burn severity should not be over-generalized to different studies performed at different geographic locations and other fire sites. The non-linear effects of topographic characteristics have been reported in several other studies (e.g., [3,66]). We considered two components of the effects of topography on burn severity, direct and indirect effects, although there are no clear-cut criteria to distinguish between these types of effects. Direct effects can be considered the physical settings, such as convex/concave terrain [75] and slope gradient [75], of a fire event. Meanwhile, indirect effects can be considered to affect the pre-fire conditions of fuel (e.g., composition, configuration, density, average stand diameter and long-term moisture) [19,76]. As Lee et al. [3] discussed, it is neither clear nor well documented in fire research why topographic variables are nonlinearly associated with burn severity. However, these two components of topographic effects are involved in the relationship between topographic characteristics and burn severity. From the perspective of direct effects, the effects of susceptible tree cover might become too great, overriding the effects of topographic characteristics [3]. Thus, burn severity could become more or less severe due to the availability of susceptible tree cover under the same topographic characteristics. Extending this rationale, additional causes of nonlinear effects of topographic variables could include the spatial variation in fuel distribution, fuel moisture, moisture content, temperature, precipitation and compounding effects [3,30,66].
From the perspective of indirect effects, topographic characteristics could alter the long-term pre-fire conditions of fuels, in turn affecting direct effects. As confirmed in many other studies, pre-fire stand conditions such as tree density, canopy coverage, average stand diameter and beetle outbreak history have notable influences on the spatial patterns of fire severity [20,26,77]. Therefore, burn severity should be understood based on the long-term interactions between topographic characteristics and fuels conditions and the spatial and temporal dimensions must be considered together; however, such a highly complex nonlinear relationship was beyond the scope of this study. Typically, the spatial distribution of fire severity is not linear with elevation and slope [78,79]. We observed similar patterns of red pine trees with elevation, slope, TWI and SRI at the study site, although we could not consider the temporal dimension. We estimated the GAM for the percentage of red pine trees with all topographic characteristics. The R 2 of the model was 0.578 and an identity link function was used to estimate the model (Table 6). Therefore, elevation, slope, TWI and SRI explained about 57.8% of the percentage of red pine trees and the effects of the topographic variables were nonlinear. Elevation (EDF = 5.82) and slope (EDF = 6.94) showed stronger nonlinear patterns, while TWI (EDF = 1.59) exhibited a weaker nonlinear pattern. These results were consistent with the findings of Miller and Urban [78], who reported nonlinear distributions Table 6. Estimated GAM. The GAM estimated about 57.8% of the variation in the percentage of red pine trees. All of the topographic variables showed significant non-linear relationships with the percentage of red pine trees. In Equation (4), the burn severity at a given location is defined by the constant, linear positive effects of red pine trees and the sum of the covariates of elevation, slope, TWI and SRI. From Table 5, the linear effects of red pine trees in Equation (4) could be defined by the constant and the sum of the covariates of topographic variables, as shown in Equation (5).

Variable
g (GAM pine ) = 58.243 + s (elevation) + s (slope) + s (TWI) + (SRI) variable (5) where g is the link function (identity) and s is the smooth function of each variable. Nonetheless, elevation showed the strongest non-linear effect on burn severity, whereas slope, TWI and SRI appeared to have weak effects on burn severity. In addition, the separate GAMs for the percentage of red pine trees with each topographic variable indicated that elevation (Adj. R 2 = 0.50) was the primary factor affecting the distribution of red pine trees, while slope (Adj. R 2 = 0.22), TWI (Adj. R 2 = 0.07) and SRI (Adj. R 2 = 0.06) were less significant factors in the distribution of red pine trees. Elevation showed the strongest nonlinear effect on burn severity, whereas slope, TWI and SRI had weak effects on burn severity. However, the role of elevation on burn severity is somewhat confusing in fire research. Some studies have reported that elevation intensifies burn severity because surface fuels at high elevations typically dry more quickly due to good drainage and increased solar exposure [19,26,80]. By contrast, other studies have reported negative effects of elevation on burn severity due to the higher productivity, vegetation density and fuel accumulation (e.g., [3,81]). However, our results suggested that the role of elevation on burn severity cannot be simplified, as suggested in previous studies and elevation can either increase or decrease burn severity at a given location depending the availability of susceptible fuels, overriding their effects and the pre-fire effects of elevation on the distribution on susceptible fuels.

Conclusions
The spatial mosaics of burn severity in burned areas have profound impacts on forest ecosystems and the response of various components of forest ecosystems to forest fires. Broadly, burn severity is determined by the forest fire triangle (i.e., fuel, topography and fire weather) and their interactions. Many studies have investigated the separate or combined effects of these factors on burns in various geographic locations. Most studies have assumed linear relationships and over-simplified the impacts of various factors on burn severity, which could mislead policy makers, managers, planners and practitioners managing fire-resilient or sustainable forests. However, the linearity of the relationships between environmental variables and burn severity has never been examined, although a few studies have reported the possible presence of nonlinear relationships between environmental variables and burn severity; therefore, the relationships between environmental variables and burn severity are not fully understood.
In this study, we examined the relationships between topographic characteristics and red pine trees with burn severity at the Samcheok fire site in Korea. We estimated LMs and GAMs between burn severity and topographic characteristics and red pine trees. After comparing the models, the nonlinear GAMs considerably outperformed the LMs in explaining burn severity. However, susceptible forest cover (i.e., Japanese red pine trees) showed linear effects on burn severity, consistent with the findings of numerous other studies, where less susceptible tree cover was associated with lower burn severity. This reinforces the importance of forest management and silviculture in lowering the availability of susceptible fuels, such as prescribed burning, thinning fuel density, planting less susceptible trees, or creating fuel breaks. Meanwhile, topographic characteristics had nonlinear effects on burn severity. In particular, elevation was the primary topographic factor affecting burn severity in a nonlinear pattern, while the nonlinear effects of slope, TWI and SRI were modest. The nonlinear effects of elevation on burn severity could be due to overriding effects of susceptible fuels over elevation effects and nonlinear effects of topographic characteristics on pre-fire fuel conditions, including the spatial distribution and availability of susceptible tree cover.
These results might only be true in moderate weather conditions, because some studies have reported that extreme weather could override the limits of fuels and topographic characteristics.
Due to a lack of fire weather data, we could not examine the interactions among fuels, topographic characteristics and weather conditions. Therefore, weather conditions should be considered when investigating the nonlinear effects of topographic characteristics and susceptible tree covers in future studies. In addition, the linear effects of red pine trees and nonlinear effects of elevation observed in this study must be validated using different tree types and geographic locations and the findings of this study should not be over-generalized for different tree cover types and geographic locations. Finally, we could not fully explain the mechanism of the nonlinear effects of elevation on burn severity. To understand the nature of the nonlinear effects of elevation, as well as other topographic variables, the same study should be replicated at several different fire sites to extract common mechanisms of the detailed effects of topographic characteristics on burn severity.