Siberia contains the only deciduous needleleaf forests (that is, larch forests, Larix spp.) underlaid by permafrost in the world, and these forests occupy 20% of the boreal forest biome (Abaimov 2010). In addition to the unique ecology of larch trees, Siberian ecosystems are distinctive because continuous permafrost occurs at subarctic latitudes, and because much of northern Siberia is underlain by thick, carbon and ice-rich permafrost deposits. Larch forests protect the underlying permafrost, thereby ensuring permafrost stability (Paulson and others 2021;Walker and others 2021;Hewitt and others 2022;Loranty and others 2024). Fires are common albeit infrequent in Siberian larch forests (fire return interval is 65-125 years (Talucci and others 2022a)), and their annual burned area is an order of magnitude greater than that of any other vegetation type in the permafrost zone (Loranty and others 2016). These fires modulate permafrost conditions and post-fire vegetation composition and forest structure (Alexander andothers 2012, 2018) and are crucial to the persistence of Siberian larch forests (Kharuk and others 2021).
Total burned area in the Siberian larch region varies considerably from year to year, with high fire activity strongly related to large-scale atmospheric conditions such as positive Arctic oscillations, which are associated with higher-than-normal temperatures in Eurasia (Balzter and others 2005;Kim and others 2020), and Arctic front jets, which bring strong winds and anomalously warm and dry conditions to northern Siberia (Scholten and others 2022). Interannual variation in burned area is also related to interannual variation in weather, including regional air temperature, soil moisture, and drought indices (Jupp and others 2006;Bartsch and others 2009;Forkel and others 2012;Ponomarev andothers 2016, 2018;Tomshin and Solovyev 2021; Descals and others 2022; Talucci and others 2022a) as well as the timing of snowmelt, which determines the period of fuel drying (Kim and others 2020;Scholten and others 2022;Talucci and others 2022a). Given that rapid increases in air temperature and earlier snowmelt have been observed across the region (Box and others 2019;Dauginis and Brown 2021), burned area has likely increased across Siberia over recent decades. However, the lack of consistent, moderate resolution satellite images over Siberia prior to 2000 combined with high interannual variability in burned area makes detecting long-term change difficult (that is, higher resolution satellites such as Landsat (30 m) can map smaller burn scars but are data-limited prior to 2000, whereas coarser resolution satellites ( ‡ 500 m) have a longer record over Siberia, but underestimate both burn area and number of fires (Talucci and others 2022a)). Nonetheless, several studies show recent increases in burned area in Siberian larch forests (Ponomarev and others 2016; Garcı ´a-La ´zaro and others 2018; Kirillina and others 2020;Tomshin and Solovyev 2021), while other analysis indicate a positive trend but no significant change in burned area (Jones and others 2022).
With continued climate change, burned area is expected to increase in Siberia (Sherstyukov and Sherstyukov 2014;Williams and others 2023) due to a concomitant increase in lightning strikes and drier fuels (Finney and others 2018;Chen and others 2021;Hessilt and others 2022). These environmental changes could lead to more ignitions (that is, number of fires), larger fires, or a combination of both. Understanding burned area changes in terms of the number of fires and fire size has significant implications for the persistence of larch forests and for permafrost stability. Large fires, for example, could increase the distance from burned areas to seed sources since larger fires tend to have a lower density of fire refugia (Talucci and others 2022b). This is important because larch seeds do not survive fire and do not form a persistent seedbank, instead depending on wind dispersal from nearby unburned trees (Abaimov 2010). Increasing distance to seed source could therefore result in lower density forests or a complete shift in functional type from trees to shrubs and/or grasses (Abaimov and Sofronov 1996;Cai and others 2013;Barrett and others 2020). Because tree density modulates soil temperature, carbon storage, and albedo in larch forests (Suzuki and Ohta 2003;Alexander and others 2012;Webb and others 2017;Kropp and others 2019;Loranty and others 2024), increasing fire size could restructure the ecosystem, with significant feedbacks to permafrost stability as well as regional and global climate.
In the short-term, an increase in burned area will reduce fuel loads and increase fuel heterogeneity across the landscape, with implications for permafrost stability (soil organic layer loss and forest edges promote permafrost thaw) (Jafarov and others 2013;Nossov and others 2013;Baltzer and others 2014;Holloway and others 2020), forest resiliency (small and isolated stands tend to be more susceptible to climate warming) (Khansaritoreh and others 2017), and future fire activity. Over the longer term, such a reduction in fuel loads/continuity could act as a negative feedback to climate change, if extreme fire weather is unable to promote fire spread due to insufficient fuel availability (Kelly and others 2013;He ´on and others 2014). At the same time, fires in Siberian larch forests are typically thought to be ignition-limited (that is, fuel is plentiful but too moist to ignite) (Kharuk and others 2021), and fuel availability does not appear to constrain large fires in southern Eurasian larch forests (Liu and others 2013;Fang and others 2015). The relative importance of fuel availability, fire weather, and other landscape controls on fire size in Siberian larch forests is unknown, thereby limiting our ability to project future fire activity.
In this study, we seek to better understand fire activity in Siberian larch forests by differentiating trends in burned area from those in fire size and the number of fires and by examining what controls each of these fire characteristics. We studied fires that occurred between 2001 and 2020 in Siberian larch forests underlaid by continuous permafrost, with the goals to understand (1) trends in burned area, fire size, and number of fires, (2) drivers of interannual variability in burned area, fire size, and number of fires, (3) drivers of spatial variability in fire size, and (4) the sensitivity of burned area, fire size, and number of fires to weather variables (precipitation and temperature) and landscape type. To study the effects of landscape type, we applied a clustering analysis to landscape variables (slope, elevation, vegetation, and water cover, and so on). Our study differs from previous studies of burned area in Siberian larch forests in that we separate the controls of and trends in fire size, number of fires, and total burned area as well as distinguish how weather variables impact burned area differently across landscape types.
We evaluated the drivers of fire characteristics in Siberian larch forests using perimeters for fires that occurred between 2001 and 2020 together with geospatial datasets of landscape, weather, and fuel characteristics. We considered the following weather variables: vapor pressure deficit, climactic water deficit, Palmer drought severity index, soil moisture, wind, and maximum air temperature in fire month, precipitation and temperature anomalies in the preceding summer, annual precipitation and temperature anomalies, and the meltwater anomaly (see Table 1 for a complete list of variables and data sources). We chose these variables, rather than metrics of fire weather (for example, the fire weather index (Field and others 2015), which are better suited to predict the probability of fire at a given location and time), to be able to disentangle the broad-scale climate processes controlling burned area, fire size, and number of fires.
To determine the drivers of interannual variability in fire characteristics, we fit linear regressions relating weather variables to annual burned area, mean fire size, and number of fires. We fit separate linear regressions to determine if the trends in these fire characteristics (annual burned area, mean fire size, and number of fires) were significant over our study period. To determine drivers of spatial variability in fire size, we fit a machine learning model relating fire size to landscape, weather, and fuel characteristics. Additionally, we applied a clustering algorithm to the landscape characteristics within fire perimeters, which identified two primary clusters: upland and lowland fires. We used the upland/lowland classification to structure subsequent analyses.
Our study region was larch forests in the continuous permafrost zone of Eurasia north of 50°N. We delineated larch forests using the 2015 European space agency climate change initiative land cover map (Defourny 2017) and delineated permafrost zones according to Obu and others (2018). The study region covers 2.9 m km 2 , which is 70% of the Eurasian continuous permafrost region.
Burned area, fire size, and number of fires were taken from an existing Landsat-derived dataset of 2001-2020 fire perimeters across Siberia (Talucci and others 2021). The fire perimeters capture a range of burn severities and likely include both high and low severity fires (that is, fires where canopy trees survive) (Talucci 2022a, b). From this dataset, we selected fires within the continuous permafrost zone (Obu and others 2018) where > 10% of the pixels (30 m) within the fire perimeter were larch dominated (Defourny 2017). The dataset we analyzed included 8886 fires across 20 years (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020).
We calculated fire size as the area within each fire perimeter, and we calculated total burned area as the sum of all fire sizes. Due to limitations of Landsat images collected in some years (that is, the Landsat 7 scanline error), it was not possible to quantify the area of unburned patches (that is, refugia) in the entire 2001-2020 fire perimeter dataset (Talucci andothers 2021, 2022b). Previous work suggests that refugia are more likely in topographic depressions and in areas with steep slopes, low tree cover, and higher elevations (Talucci and others 2022b), which suggests that our overestimation of burned area could be higher in uplands than in lowlands.
To understand the relative contribution of different fire sizes to total burned area, we divided these fires into four size classes: small (0-10 K ha), medium (10-100 K ha), large (100 K-1 M ha), and mega > 1 M ha). For each year (2001-2020), we quantified the number of fires and the burned area in each size class.
To determine if the burned area, mean fire size, or number of fires changed over the 20-year study period, we fit linear models relating each of these variables to fire year. We fit separate models for uplands, lowlands, and the entire study region. We used the Pearson correlation coefficient to assess the association between annual burned area and number of fires and mean fire size.
We considered multiple weather variables (in separate, simple linear regressions; see below) as A subset of these variables (weather conditions) were used in the analysis of interannual variability in mean fire size, number of fires, and burned area (see Methods for more details).
potential drivers of the interannual variation in burned area, mean fire size, and number of fires.
When modeling burned area and the number of fires, each weather variable was averaged across the entire study region in each year, because we predicted that the burned area and number of fires in the study region would depend on regional weather. In contrast, when modeling mean fire size, we first averaged weather variables across the pixels within each fire perimeter and then averaged the weather variables across all fire perimeters in a given year, because we predicted that fire size would depend on the weather conditions at the fire locations (rather than the regional mean weather).
The same weather variables were used for all three response variables (burned area, mean fire size, and number of fire), but the weather variables were calculated in different ways with respect to spatial averaging (as explained above) and timing (see below).
Weather variables included in our analysis were: vapor pressure deficit, wind speed, climatic water deficit, Palmer Drought Severity Index, monthly maximum air temperature, soil moisture, annual precipitation, annual temperature, and precipitation and temperature in the preceding summer. For models of burned area and number of fires, vapor pressure deficit, wind speed, climatic water deficit, Palmer Drought Severity Index, monthly maximum temperature, and soil moisture were averaged across the months of May-August ( 85% of fires occur between these months) in each year. For the mean fire size model, these variables were evaluated during the month of the first day of the fire. Annual precipitation and temperature variables were calculated as the sum of precipitation or temperature from snow-off (that is, the last snowon date) in the previous year to snow-off in the fire year. Precipitation and temperature in the preceding summer were calculated as the sum of the precipitation or temperature from snow-off in the previous year to the first snow-on in the previous year. Annual and preceding summer temperature and precipitation data came from the ERA5-Land hourly dataset ( 9 km) (Mun ˜oz 2019), snow cover timing was derived from the daily MODIS snow cover product ( 500 m) (Hall and Riggs 2016), and other weather variables, including monthly maximum air temperature, were obtained from the monthly TerraClimate dataset ( 4 km) (Abatzoglou and others 2018). We included variables from both the ERA5-Land dataset and the TerraClimate dataset because while the ERA5-Land dataset has higher temporal resolution (hourly vs. monthly) and is therefore more suitable to char-acterize fire weather, it does not include integrated metrics of temperature and precipitation (for example, vapor pressure deficit, climactic water deficit, Palmer drought severity index), which were also important in our analysis.
We fit simple linear regression models that related annual burned area, annual mean fire size, and annual number of fires to each weather variable. We chose this approach because collinearity among predictor variables made model selection unreliable using a multiple regression approach, and because the size of the dataset (n = 20 years) was too small to employ more complicated methods without overfitting the models.
Prior to all linear regression and correlation analyses described above, burned area was square root transformed and the number of fires and mean fire size were log transformed to help meet the assumptions of linearity, normally distributed errors, and homoscedasticity. We visually assessed the scatterplots of the correlated variables and found no obvious departures from these assumptions. We implemented the Shapiro-Wilk normality and Breusch-Pagan tests to test the residuals of our linear models for normality and heteroskedasticity, respectively. The sample size for these analyses was n = 20 years, which we judged to be insufficient for more complex analyses (for example, nonlinear models with non-normal errors). We implemented the correlation and regression analyses using the cor, cor.test, and lm functions, respectively, in R (R Core Development Team 2023). We report results for all regressions where the model p-value was less than 0.05.
We explored different weather, fuel, and site variables as potential drivers of fire size, which are summarized in Table 1 and described below. The original pixel resolution varied depending on the original dataset (see Table 1), but all datasets were resampled to 30 m using the nearest neighbor method. For most variables, pixel-level values were averaged across all pixels within each fire perimeter. The exceptions to this were the proportion of pixels classified as larch or water within each fire perimeter and categorical variables (for example, Arctic/subarctic), which were calculated (proportions) or obtained (categories) for each fire perimeter. All geospatial processing was conducted in Google Earth Engine (Gorelick and others 2017).
Weather variables included in our fire size analysis were either evaluated in the month of the fire (climatic water deficit, Palmer drought severity index, vapor pressure deficit, windspeed, and monthly maximum temperature from TerraClimate ( 4 km) (Abatzoglou and others 2018)), or as the anomaly from the 20-year (2001-2020) mean (meltwater, annual temperature, annual precipitation, summer temperature, and summer precipitation, calculated from the ERA5-Land hourly dataset ( 9 km)). To derive weather anomalies, we first calculated the total amount of water in snowmelt and the sums of annual precipitation, summer precipitation, annual temperature, and summer temperature in each year (2001-2020) at each pixel ( 9 km; native resolution of the dataset) across the study region. Anomalies were computed as the difference between these variables at the location of the fire in the fire year and the twenty-year mean at the location of each fire. To derive annual meltwater, we summed the amount of water in snowmelt from January through July. To derive annual precipitation and temperature, we summed precipitation or temperature from snowoff in the previous year to snow-off in the current year. To derive precipitation and temperature in the previous summer, we summed precipitation or temperature from snow-off in the previous year to snow-on in the previous year. Snow-on and snowoff dates were derived from the daily MODIS snow cover product ( 500 m) (Hall and Riggs. 2016).
Site characteristics included in our fire size analysis were long-term mean annual precipitation and mean annual temperature (Hijmans and others 2005), latitude, percentage of sand within the soil (Poggio and others 2021), terrestrial ecozone (for example, mountain tundra, taiga) (Olson and others 2001), whether the fire occurred in the Arctic or subarctic (Talucci and others 2021), the proportion of water pixels within a fire, and topographic variables derived from digital elevation models (slope, ruggedness relative to the surrounding 7 ha landscape, and elevation). Water pixels were defined as any pixel identified by Pickens and others (2020) as containing water (that is, experiences seasonal inundation, is permanent water, or is wet with high frequency) over the Pickens and others (2020) dataset period . Topographic variables were derived from the Arctic digital elevation model (DEM; 2 m) (Porter and others 2018) for pixels above 60°N and the NASA DEM ( 30 m) (NASA JPL 2020) for pixels below 60°N (the Arctic DEM does not extend below 60°N and the NASA DEM does not extend above 60°N). We included the Arctic/subarctic category because previous work indicated a different fire regime between the two regions, potentially related to permafrost depth and forest stand structure (Talucci and others 2022a).
We considered multiple above-and belowground fuel characteristics to account for the fact that fires can be surface fires (that is, fueled by aboveground ground layer vegetation), ground fires (fueled by soil organic matter), and/or crown fires (fueled by tree crowns), burning both aboveand belowground vegetation (Webb and others 2024). Specifically, we included the proportion of larch pixels within a fire (larch pixels were identified from Defourny (2017)), soil carbon density in the top 5 cm of the soil profile (a proxy for soil organic matter fuel loads) (Poggio and others 2021), tree canopy cover in 2000 (a proxy for aboveground fuel loads (Alexander and others 2024)) (Hansen and others 2013), tree cover connectivity, soil moisture in the month of fire (Abatzoglou and others 2018), and the first day of the year associated with the fire ('first burn day') (Talucci and others 2021). We defined tree cover connectivity as the proportion of the area within the fire that is both connected forest (defined as spatially contiguous groups of pixels connected at one or more edges) and where forest cover is > 30% (Hansen and others 2013). We considered the proportion of larch separately from tree cover because larch species have traits associated with fire resistance such as high leaf moisture, thick bark, and self-pruning of lower branches that control fire intensity (and therefore also likely control fire spread) differently from other co-occurring tree species (Wirth 2005;Rogers and others 2015). We considered the first burn day as a fuel characteristic because as the active layer progressively thaws over the course of the summer, more organic soil is available for combustion (Turetsky and others 2011).
We fit a histogram-based gradient boosting regression tree (HGBRT) model (a machine learning algorithm) that related log transformed fire size to site characteristics, fuel characteristics, and weather variables (see above and Table 1). We chose the HGBRT approach because regression trees are an easily interpretable method of determining variable importance and ensemble methods such as boosting generally produce better models (lower bias and variance) than single tree methods (Elith and others 2008). We fit the HGBRT using the Python-based Scikit-learn library (Pedregosa and others 2011), and optimal model hyperparameters were determined using grid search and tenfold cross-validation (Kohavi 1995;Elith and others 2008). We selected the best performing model and evaluated its performance using tenfold cross-validation repeated 100 times; the mean model R 2 (standard deviation) was 0.46 (0.03). This final model was used for subsequent analyses of feature importance and partial dependence of features.
Using the HGBRT model, we employed permutation importance (Pedregosa and others 2011) to quantify the relative importance of each potential driver of spatial variability in fire size. Permutation importance randomly shuffles the value of each explanatory variable among all fires and determines the resulting drop in average R 2 value. To quantify the relative importance of each group of explanatory variables (that is, site characteristics, fuel characteristics, and weather conditions), we permutated each group of explanatory variables together, rather than individually. Individual and grouped variable permutation importances were repeated 100 times and were implemented using the Scikit-learn (Pedregosa and others 2011) and rfpimp (Parr and Turgutlu 2018) libraries, respectively.
To understand the shape and direction of the relationships between fire size and explanatory variables, we used the HGBRT model to derive the marginal effect of the two most important variables in each group (site characteristics: proportion of water pixels in the fire perimeter and slope; fuel characteristics: proportion of larch pixels in the fire perimeter and first burn day; weather conditions: climactic water deficit and vapor pressure deficit) on fire size. Marginal effects are calculated by generating a fire size prediction for each observation in the dataset while varying the value of the intended variable and keeping all other explanatory variables as is. The marginal effect of a value of the intended variable is the average fire size prediction at that value, and the range of marginal effects is plotted in a partial dependence plot (Hastie and others 2009).
Initial analyses revealed a sharp directional change in marginal effect between 0 and 0.004 for the proportion of water within the fire perimeter and between 0.996 and 1 for the proportion of larch pixels within the fire perimeter. Because it is unlikely that such small changes in water or larch presence could, by themselves, have strong effects on fire behavior, these apparent effects likely reflect the presence of unmeasured variables (or errors in the measured variables) rather than physical processes. For example, it seems unlikely that the difference between the water fractions 0 and 0.004 would, by itself, cause a significant change in fire behavior, whereas it seems plausible that this apparent effect in reality reflects a meaningful dif-ference in site conditions (for example, thin rocky soils vs. deeper soils) that was not captured by the available data (Table 1). Given the ambiguous interpretation of the marginal effects over these small intervals (0-0.004 for proportion of water; 0.996-1 for proportion of larch), we have excluded these intervals from the partial dependence plots shown in the main text; for completeness, we show the non-truncated partial dependence plots in Figure S1.
Grouping Data into Wet/Dry Years, Cool/ Hot Years, and Landscape Types
We also studied the sensitivity of burned area, fire size, and number of fires to categorical groupings of weather (low/high precipitation/temperature years) and landscape variables (upland/lowland ecosystems). To group weather variables, we averaged the annual air temperature across the fire perimeters of each year and then classified years as low (below the median) or high (above the median) temperature. We similarly classified years as wet (above median) or dry (below median) based on annual precipitation (averaged across fire perimeters of each year).
To group fires into different landscape types, we applied the k-means clustering algorithm to elevation, ruggedness, slope, tree cover, soil carbon density, and percent sand within the soil (see above and Table 1 for variable definitions and data sources). These variables were selected because they were important drivers of fire size (see Results) and because they are defined for individual pixels, which simplified applying the clustering algorithm at the regional scale (see below). We did not include variables that could only be calculated for multiple pixels (for example, proportion of water pixels within a fire) because it was not straightforward to include such variables in the regional analysis described below. We implemented the kmeans clustering algorithm using the KMEANS function in the fdm2id library (Blansche ´2023) in R (R Core Development Team 2023). The optimal number of clusters, based on the silhouette method, was two. Based on the mean values of the characteristics of the two clusters (Table 2), we named the clusters 'upland' (higher elevation, steeper slopes, lower prevalence of water, and lower soil carbon density) and 'lowland' (lower elevation, gentler slopes, higher soil carbon density, and higher prevalence of water). The terms 'upland' and 'lowland' are widely used to describe landscape position in the boreal forest biome (for example, Chapin and others 2010; Eichhorn 2010; Jorgenson and others 2022), but these terms may be ambiguous in some cases. For example, depending on the values of the six classification variables (elevation, ruggedness, slope, tree cover, soil carbon density, and percent sand), some fires on plateaus may be classified as upland, whereas other fires plateaus may be classified as lowland. Thus, we use the terms 'upland' and 'lowland' for convenience, but these terms may not be suitable for every pixel in our study and may not have a simple correspondence to other topographic terms.
To calculate the regional fraction of upland and lowland ecosystems, we applied the same k-means clustering algorithm described above to the entire study region. We first extracted the variables necessary for the cluster analysis (elevation, ruggedness, slope, tree cover, soil carbon density, and percent sand within the soil) from 100,000 randomly selected 30 m pixels from across the study region. (We sampled a subset, rather than the entire study region, to reduce the computational demands of our analysis.) We then used the k-means clustering results based on fires to classify each of these 100,000 pixels as either upland or lowland, which provided a region-wide estimate for the fractions of upland and lowland pixels.
There was no statistically significant trend in burned area, mean fire size, or number of fires in uplands, lowlands, or across the entire study area (p > 0.05 for all trends; Figure 1). Interannual variation in burned area was large, with the largest Siberian Larch Forests fire year burning more than 36 times as much area as the smallest fire year. Mean fire size also varied substantially from year to year; the largest annual mean fire size was 6.5 times greater than the smallest. Similarly, the year with the most fires had nearly 7 times more fires than the year with the fewest. Annual burned area was strongly related to both the number of fires (r = 0.86; p < 0.01) and the mean fire size (r = 0.87; p < 0.01) (Figure 2), although the number of fires and the mean fire size were considerably less correlated (r = 0.54; p = 0.01). The relative importance of different fire sizes varied across years, with fire sizes over 1 M ha occurring only in 2002 (Figure 3). When all years were combined, small fires (< 10 K ha) accounted for 88% of all fires (n = 7,836) and 25% of the area burned, medium fires (10-100 K ha) accounted for 11% all fires (n = 938) and 37% of the area burned, large fires (100 K-1 M ha) accounted for 1% of all fires (n = 110) and 35% of the area burned, and mega fires (> 1 M ha) accounted for less than 1% of the number of fires (n = 2) and 4% of the area burned. Interannual variation in burned area and number of fires were most strongly linked to regional vapor pressure deficit, but precipitation, temperature, and integrated precipitation/temperature metrics (that is, Palmer Drought Severity Index, climatic water deficit, and soil moisture) were also important predictors of burned area and number of fires (Table 3). The only statistically significant predictor of interannual variation in mean fire size was precipitation in the preceding year (Table 3).
Site characteristics were the most important predictors of fire size, followed by fuel characteristics and weather conditions (Figure 4). Two of the three most important predictors of fire size (proportion of water within a fire perimeter and slope; Figure 4) were also important in distinguishing upland from lowland ecosystems (Table 2). Fires tended to be larger in areas with a small amount of water and with shallow slopes (Figure 5). The proportion of larch pixels within a fire and the first day of the fire ('first burn day') was the most important fuel characteristics, with higher proportions of larch pixels and earlier first burn days associated with larger fires (Figure 5). Climatic water deficit and vapor pressure deficit were the most important weather variables, with higher deficits leading to larger fires (Figure 5).
Upland ecosystems (that is, those with steeper slopes, rugged terrain, and higher elevation) occupied 31% of the study region, but contained 44% of the fires, suggesting that uplands are more likely to ignite than lowlands (Table 2). However, uplands accounted for a disproportionately small percentage of the burned area (20%) because, on average, fires in uplands were 1/3 the size of fires in lowlands. Lowland ecosystems (that is, those with gentle slopes, higher proportion of water, and higher soil carbon density) occupied the other 69% of the study region and supported a comparatively small number of fires (56%) but sustained a disproportionately large percentage of the burned area (80%) (Table 2).
In uplands, fire characteristics were more sensitive to temperature than to precipitation, with high-temperature years associated with a 35% increase in mean fire size, 45% increase in the number of fires, and 95% increase in burned area (Figure 6). Greater burned area in high-temperature years was due to increases across a wide range of fire sizes (Figure 7). In comparison, low-precipitation years were associated with smaller effects on mean upland fire characteristics (Figure 6), although the largest upland fires occurred in lowprecipitation years (Figure 7).
In contrast to uplands, lowlands were more sensitive to precipitation than to temperature, with low-precipitation in the preceding year associated with a 65% increase in mean fire size, 79% increase in number of fires, and 194% increase in burned area (Figure 6). Temperature effects were weaker than precipitation effects in lowlands but were still substantial for burned area (57%) and mean fire size (50%). In years where the previous year had low precipitation, there were increases in lowland burned area across all fire sizes (Figure 7). In contrast, greater lowland burned area in years where the previous year had high temperatures was primarily due to increases in the larger fire sizes (Figure 7). The largest two fires (that is, mega fires; > 1 M ha) occurred in years with both low Estimates are slopes from simple linear regressions where the response and explanatory variables were both standardized to unit variance. To meet the assumptions of normality and homogeneity of variance, burned area was square root transformed and mean fire size, number of fires, soil moisture, and summer temperature were log transformed prior to standardization. Only significant regressions (p < 0.05) are reported. The sample size for all regressions is 20 years (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). See Table 1 for explanation and data sources of variables. *Averaged across the months of May-August for each year.
preceding year precipitation and high preceding year temperature (Figure 7). The entire study region (uplands and lowlands combined) was more sensitive to precipitation than to temperature, with low precipitation in the preceding year associated with 76% larger fires, 31% more fires, and 131% more burned area when precipitation was high in the preceding year (Figure 6). These trends in fire size and burned area are largely driven by the sensitivity of lowlands to precipitation (Figs. 6 and 7), since lowland fires account for the largest fires and 80% of the total burned area (Table 2).
Spatial variation in fire size was most strongly related to landscape position. Whereas there were more fires in landscapes with steeper slopes, greater terrain ruggedness, and higher elevation (that is, uplands), fires were larger in flatter, lower elevation areas with lower terrain ruggedness (that is, lowlands). These results likely reflect different fuel characteristics in uplands and lowlands. Uplands tend to be better drained, leading to drier surface fuels that more readily ignite, but fuel is also more discontinuous with more natural firebreaks (for example, rocky outcrops, creeks), limiting fire spread (Sofronov and Volokitina 2010). A larger number of fires in uplands could also reflect the fact that steep slopes and rugged terrain promote spot fires (Storey and others 2020), with fire spotting resulting in additional, distinct fire perimeters beyond the primary fire perimeter.
Fires in Siberian larch forests are typically surface fires, fueled by the moss and lichen matrix on the forest floor and the deep undecomposed soil organic layer rather than the forest canopy (Sofronov and others 2000;Kharuk and others 2021). Because moss thrives in wet environments, lowlands tend to have higher fuel loads with fewer firebreaks (Sofronov and others 2000), which means that, once ignited, these fires can burn over larger areas, particularly following low-precipitation years when the fuel is relatively dry. In well-drained areas like uplands, water has a shorter residence time, so surface fuels are unlikely to maintain high water content for prolonged periods (for example, multiple weeks or months), even during high-precipitation years. The rapid draining of uplands may explain why upland fire size is mostly insensitive to precipitation.
In lowlands, our results suggest that wet conditions (that is, high precipitation in the preceding year) protect surface fuel from excessive burning, but that dry conditions may lead to increases in fire size across all fire size classes. This corroborates earlier work that demonstrated that soil moisture tends to act as an 'on/off switch,' with large burned areas not possible if soil moisture is moderately high (Bartsch and others 2009). Low precipitation may dry out natural firebreaks such as typically wet mossy bogs, streams, and shallow rivers, allowing fire to spread farther in dry years (Sofronov and Volokitina 2010). At the same time, low precipitation also dries out forest floor fuels, increasing the flammability surface fuels. Our results also corroborate Forkel and others (2012), who showed that high burned area is related to low soil moisture conditions in the previous year. This lagged effect is likely due to two co-occurring mechanisms. First, because the ground is frozen during snow melt, meltwater does not infiltrate the soil profile, so fall soil moisture conditions largely control spring soil moisture conditions (Sofronov and others 2000). Second, because water has a high specific heat, the previous year soil moisture affects the rate of soil thaw in the spring (Sofronov and others 2000), with low moisture years thawing organic layer fuels earlier than high moisture years.
Across the entire study region (uplands and lowlands combined), fires that were initiated in the spring tended to be larger, with fire size decreasing over the course of the season. This is the opposite pattern of observations of fire size in North American boreal forests, where later season fires are larger because more of the active layer is thawed, making more ground fuels available for combustion (Turetsky and others 2011). However, our results are consistent with field-based observations in Siberian larch forests, where early season 'runaway' fires fueled mostly by the litter layer (since the soil organic layer is not thawed) are common (Sofronov and Volokitina 2010;Kharuk and others 2021).
Most of the precipitation in Siberia occurs during the summer, particularly in July and August (Kostrova and others 2020;Han and Menzel 2022). This mid-to late-season precipitation could reduce fire size by expanding the presence of natural fire breaks (for example, wet mossy bogs, streams), and by moistening lichens and mosses (Mallen-Cooper and others 2021) and the underlying soil organic layer that make up the majority of the fuel load.
Larch forests have the largest relative burned area (per-unit land area) of any forest type in Siberia (Kharuk and others 2021). Similarly, we found that fires with a higher proportion of larch pixels tended to be larger, likely due to species traits and stand characteristics that promote fire spread. The low canopy closure characteristic of larch forests, for example, allows wind to freely penetrate the canopy, advancing the fire front (Sofronov and Volokitina 2010). Because larch trees drop their needles each fall, the presence of larch trees creates a low bulk density fuel bed of larch needles, which, combined with the underlying moss layer, promotes flammability. Additionally, larch presence directly affects understory species composition, with implications for both fuel loads and fuel moisture (Loranty and others 2018;Paulson and others 2021). In particular, larch presence increases moss abundance, a key fuel source when alive and dead (Alexander and others 2020;Paulson and others 2021). An increasing proportion of larch pixels within a fire perimeter may also reflect higher ground fuel connectivity, which may be somewhat independent of our canopy-derived metric of fuel connectivity.
Ultimately, fire size is influenced by the interaction of multiple factors, many of which do not vary on human timescales (for example, topography) or are stochastic processes (for example, wind speed, lighting strikes). Projecting future fire regimes requires identifying the climate signal within the noise of these other factors. Of our studied variables, precipitation in the preceding year was the only climate variable that could explain interannual variation in mean fire size across the study region, with lower precipitation leading to larger fires. While temperature is an important driver in uplands, precipitation is more important in lowlands. Because lowlands account for more fires and have a larger mean fire size than uplands, interannual variability in burned area across the entire region is largely determined by lowlands.
Fires in northern Siberian larch forests are primarily lightning ignited, while at the southern extent of the forest where population density is higher, anthropogenic ignitions are more common (Kirillina and others 2020;Kharuk and others 2021;Xu and others 2022). The number of fires in any year is thus a function of the number of lightning strikes, anthropogenic activity, and the susceptibility of fuel to ignition sources. We found that interannual variability in the number of fires was primarily driven by drought indices such as vapor pressure deficit, which impacts both the number of lightning strikes and fuel flammability (Sedano and Randerson 2014;Scholten and others 2022). Specifically, hot and dry conditions (that is, high vapor pressure deficit) are associated with ignition and fire spread efficiency (Sedano and Randerson 2014;Hessilt and others 2022). While we did not account for human activity in our models, fuel conditions are agnostic to the ignition source, and the hot and dry conditions that amplify fire spread in lightning-ignited fires would also increase ignition and fire spread probabilities in human-ignited fires.
Correlations between interannual variability in region-wide burned area and environmental variables were strongest for vapor pressure deficit and other integrated metrics of moisture and temperature, consistent with previous analyses of burned area in Siberia (Balzter and others 2005;Ponomarev andothers 2016, 2018;Talucci and others 2022a). Annual burned area was strongly and positively correlated with both the annual mean fire size and the annual number of fires. However, while the highest burned area occurred in years with both a high number of fires and a large mean fire size, mean fire size showed only a moderate correlation with the number of fires. This may reflect the randomness of where lightning strikes occur combined with the importance of ignition location to eventual fire size (for example, high lightning years may not result in large burned areas if the majority of ignitions occur in uplands). Additionally, interannual variability in the number of fires and fire size were best predicted by different environmental variables (current year fire vapor pressure deficit and previous year precipitation, respectively), making it less likely that optimal conditions to support both large fires and a large number of fires co-occur.
A growing body of evidence, including this study, demonstrates that burned area in Siberian larch forests is highly dependent on temperature, precipitation, and the combination of the two (Jupp and others 2006;Bartsch and others 2009;Forkel and others 2012;Ponomarev andothers 2016, 2018;Tomshin and Solovyev 2021;Descals and others 2022;Scholten and others 2022;Talucci and others 2022a). While models project both warming air temperatures and increasing precipitation across Siberia, future precipitation projections are more uncertain than temperature projections (Van Der Wiel and Bintanja 2021). For example, some climate models project an increase in winter and fall precipitation with little to no change in the summer (Cai and others 2024), but observations across Siberia show that most of the increase in precipitation over the past 70 years occurred during the summer (Wang and others 2021). Accurately understanding the timing of precipitation will be important for projecting future fire size. Increased snowfall, for example, might not directly impact the fire regime if snowmelt runs off and does not contribute to soil moisture conditions (Sofronov and others 2000). On the other hand, increasing summer precipitation could impact the fire regime because summer precipitation directly affects ground fuel moisture conditions and therefore fire behavior. Similarly, there is considerable uncertainty in projections of future soil moisture conditions, with some models projecting wetting from increased precipitation and others projecting drying from increased evapotranspiration, and all models lacking key processes that determine soil moisture conditions in permafrost systems (Andresen and others 2020).
Given that warming temperatures lead to Arctic wetting (Box and others 2019;McCrystall and others 2021), our results suggest that climate change may have opposite effects on the number of fires and fire size. Specifically, the number of fires increased with vapor pressure deficit and other metrics of drought, which are expected to increase with climate change (Yuan and others 2019), suggesting that Siberian larch forests could experience more fires in the coming century. On the other hand, mean fire size was negatively related to precipitation, and the expected increases in precipitation could result in a decrease in mean fire size. However, interannual variability in precipitation may increase, with some very dry years, even as mean annual precipitation increases (Pendergrass and others 2017). This precipitation variability could lead to extreme fire sizes. Climate change over future decades could therefore lead to more fires, smaller fires on average, and more variable fire sizes in Siberian larch forests. Ultimately, future trends and interannual variability in total burned area will be determined by multiple factors, including the degree of warming, the magnitude and seasonal timing of precipitation and temperature change, and interannual climate variability.