Geographically divergent trends in snow disappearance timing and fire ignitions across boreal North America

The snow cover extent across the Northern Hemisphere has diminished, while the number of lightning ignitions and amount of burned area have increased over the last 5 decades with accelerated warming. However, the effects of earlier snow disappearance on fire are largely unknown. Here, we assessed the influence of snow disappearance timing on fire ignitions across 16 ecoregions of boreal North America. We found spatially divergent trends in earlier (later) snow disappearance, which led to an increasing (decreasing) number of ignitions for the northwestern (southeastern) ecoregions between 1980 and 2019. Similar northwest–southeast divergent trends were observed in the changing length of the snow-free season and correspondingly the fire season length. We observed increases (decreases) over northwestern (southeastern) boreal North America which coincided with a continental dipole in air temperature changes between 2001 and 2019. Earlier snow disappearance induced earlier ignitions of between 0.22 and 1.43 d earlier per day of earlier snow disappearance in all ecoregions between 2001 and 2019. Early-season ignitions (defined by the 20 % earliest fire ignitions per year) developed into significantly larger fires in 8 out of 16 ecoregions, being on average 77 % larger across the whole domain. Using a piecewise structural equation model, we found that earlier snow disappearance is a good direct proxy for earlier ignitions but may also result in a cascade of effects from earlier desiccation of fuels and favorable weather conditions that lead to earlier ignitions. This indicates that snow disappearance timing is an important trigger of land–atmosphere dynamics. Future warming and consequent changes in snow disappearance timing may contribute to further increases in western boreal fires, while it remains unclear how the number and timing of fire ignitions in eastern boreal North America may change with climate change.

1 Introduction

Snow cover across boreal and Arctic ecosystems is an important driver of regional hydrological cycles and the global energy balance (Swenson and Lawrence, 2012; Li et al., 2017). With climate warming, spring snow cover has decreased 11 % per decade over the Northern Hemisphere since the 1970s (Déry and Brown, 2007; Brown and Robinson, 2011). Changes in snow cover and sea ice led to a substantial decrease in the cryosphere radiative forcing across the Northern Hemisphere of around 0.5 W m⁻² from 1979 to 2008, which warms the regional and global climate (Flanner et al., 2011; Groisman, et al., 1994). This feedback contributes to accelerated warming in the northern high latitudes (Anisimov et al., 2007; Rantanen et al., 2022); however, changes in snow cover are heterogeneous across the Northern Hemisphere (Bormann et al., 2018; Suzuki et al., 2020). Over boreal North America, changes in snow cover timing show a long-term spatial divergence between earlier (later) snow disappearance timing over western (eastern) boreal North America between 1972 and 2017 (Chen et al., 2016; Bormann et al., 2018). The divergent changes in snow cover will likely have important impacts on ecosystem functioning in boreal forest and Arctic tundra (Post et al., 2009; Buermann et al., 2013) and may be the result of persistent changes in atmospheric circulations (Jain and Flannigan, 2021).

Simultaneously, over the last 2 decades, large parts of western boreal North America have experienced a rise in the number of lightning fire ignitions and burned area (Hanes et al., 2019), driven by increases in dry fuel availability (Abatzoglou et al., 2016; Hessilt et al., 2022), favorable fire weather (Sedano and Randerson, 2014), and an increase in the number of lightning strikes (Veraverbeke et al., 2017). Fire is the most widespread ecosystem disturbance in boreal North America, and these increasing trends in fire occurrence are expected to continue in the future (Flannigan et al., 2005; Balshi et al., 2009; Chen et al., 2021; Phillips et al., 2022). Early snow disappearance has previously been linked to large fires in the western United States as a consequence of longer periods of fuel drying (Westerling et al., 2006). Dry fuel availability is a prerequisite for fire ignitions (Abatzoglou et al., 2016; Hessilt et al., 2022) and may further enable rapid fire growth, thereby resulting in larger fires (Sedano and Randerson, 2014; Veraverbeke et al., 2017). The relationships between snow disappearance timing and fire behavior characteristics, such as fire ignitions and size, may vary across boreal North America and remain poorly understood (Hanes et al., 2019).

In recent years, early snow disappearance after warm winters has been linked to summer heatwaves and severe fire seasons over Siberia (Gloege et al., 2022; Scholten et al., 2022). Warm winter extremes can substantially impact ecosystem functioning until deep into the subsequent growing season (Zona et al., 2022). Early snow disappearance induces an early vegetation green-up because of early peaks in soil moisture (Gloege et al., 2022) but decreased late-season vegetation productivity (Buermann et al., 2013; Miles and Esau, 2016; Graham et al., 2017). The enhanced evapotranspiration can lead to soil desiccation and result in increased sensible heat flux later in spring (Gloege et al., 2022). This enhances atmospheric warming and drying through limited evaporative cooling (Seneviratne et al., 2010). In turn, positive geopotential height anomalies and a persistent atmospheric ridge can form (Cohen et al., 2014; Tang et al., 2014) and promote atmospheric blocking events that create favorable weather conditions for fire ignition and spread (Coumou et al., 2018; Jain and Flannigan, 2021; Scholten et al., 2022). Simultaneously, destabilization of the atmosphere increases the occurrence of convective thunderstorms and lightning (Chen et al., 2021). The increases in cloud-to-ground lightning strikes can potentially increase the likelihood of igniting dry fuels (Hessilt et al., 2022). Nonetheless, the influence of a divergent snow cover trend across boreal North America on weather, fuel dryness, and ignition timing has previously not been studied and may exhibit divergent responses to changes in the snow cover.

Earlier snow disappearance may also lead to an earlier start of the fire season, thereby lengthening and intensifying the boreal fire season (Flannigan et al., 2005; Bartsch et al., 2009; Veraverbeke et al., 2017). Defining the fire season length is not straightforward and different methods have been used to quantify the length of the boreal fire season. The fire season length has been estimated using fire weather indices as proxies for fire activity (Wotton and Flannigan, 1993; Flannigan et al., 2016). Other studies have estimated the fire season length using long-term government records, which are prone to temporal changes in accuracy and uncertainties (Hanes et al., 2019). Daily fire monitoring using the polar-orbiting Moderate Resolution Imaging Spectroradiometer (MODIS) sensors allows accurate definition of the fire season based on observed fire activity since the 2000s (Justice et al., 2002; Giglio et al., 2016, 2018). Given that the MODIS record dates back until the early 2000s, it may be possible to infer changes in fire season length across boreal North America during this period.

Here, we investigated relationships between snow disappearance and early-season ignition timing across boreal North America between 2001 and 2019. In addition, we evaluated the influence of ignition timing on fire size and assessed temporal changes in snow disappearance timing and the number of ignitions since 1980. Through satellite-derived estimates, we derived the length of the snow-free and fire seasons and assessed the influence of the length of the snow-free season on fire season length. Early ignition timing was modeled as a function of snow disappearance timing and meteorological and fire weather conditions using a linear mixed-effect model to investigate potential cascading effect of earlier snow disappearance timing. Finally, we assessed the interactions between snow disappearance timing and meteorological and fire weather conditions when modeling ignition timing through a piecewise structural equation model.

Figure 1 (a) Trend in snow disappearance timing between 2001 and 2019 derived from Moderate Resolution Imaging Spectroradiometer (MODIS) for the study domain overlaid by second-level ecoregions (US EPA, 2015) and (b) the mean annual ignition density per 100×100 km grid cells between 2001 and 2019. All pixels exceeding the average pixel snow disappearance timing ±3 standard deviations were excluded and set to not applicable (NA: gray).

2 Methodology

2.1 Study domain

The study domain includes Alaska, USA, and the majority of Canada (9.17×10⁶ km²) excluding the Canadian Arctic Archipelago and is divided into 16 ecoregions (Omernik, 1987, 1995) (Fig. 1). We used the second-level ecoregions for subcontinental comparisons (McCoy and Neumark-Gaudet, 2022). We included 14 ecoregions but further divided the Softwood Shield and Taiga Shield into eastern and western ecoregions due to their large longitudinal gradients, resulting in 16 different ecoregions in our study (Fig. 1 and Table S1 in the Supplement). The Softwood Shield was divided in accordance with the third-level ecoregion division and the Taiga Shield was split into two subregions east and west of Hudson Bay (Baltzer et al., 2021) (Fig. 1). The northernmost ecoregions (the Arctic Cordillera, Northern Arctic, and Southern Arctic) were excluded as they included very few ignitions. The southern parts of the Cold Deserts, Marine West Coast Forest, Mixed Wood Shield, and Western Cordillera were cropped out as they were not covered by the Arctic–Boreal Vulnerability Experiment Fire Emission Database (ABoVE-FED; Potter et al., 2023) (Fig. 1). Our study domain thus included Arctic tundra, boreal forest, and temperate ecosystems between northwestern Alaska and southeastern Canada, hereafter referred to as “boreal North America”.

2.2 Snow disappearance timing

We retrieved snow disappearance timing at 463 m resolution from the MODIS daily composite snow cover product MOD10A1 collection 6 between 2001 and 2019 (Hall and Riggs, 2016). This product computes the normalized difference snow index (NDSI) ranging from −1 to 1 from visible and shortwave infrared spectral data. The relationship between NDSI and estimated fractional snow cover from higher-resolution snow cover data from the Landsat Enhanced Thematic Mapper-plus (30 m) has previously been proven robust over large areas such as boreal North America (Salomonson and Appel, 2004). This allowed us to use NDSI as a proxy for fractional snow cover. We identified the Julian calendar day of snow disappearance timing as the first day a pixel had less than or equal to 15 % snow cover for a minimum of 14 consecutive days (Verbyla, 2017). We also tested a threshold of snow cover less than or equal to 15 % for a minimum of 7 consecutive days but found little difference between the two thresholds (Fig. S1, Table S2). Pixels that had burned, contained persistent cloud cover, water, or perennial snow cover (more than 250 d yr⁻¹), or contained less than or equal to 15 % snow cover for fewer than 14 consecutive days were excluded from the analysis. Pixels with values exceeding a pixel-specific threshold (average snow disappearance timing in 2001–2019 ±3 standard deviations) were regarded as outliers and excluded from the analysis. The snow disappearance timing was determined between 1 February and 31 July. We opted for a large potential range in snow disappearance timing because of the large latitudinal and thus climatological range present in the study domain (Fig. 1). To retrieve the first day of snow cover, we used the reversed method where the first day on which at least 15 % of the pixel was snow-covered for a minimum of 14 consecutive days was set to the first day of snow cover. This was determined between 1 August and 31 December. We modified the code from Armstrong et al. (2023) to compute the snow disappearance timing in Google Earth Engine.

In complement to the MODIS snow cover product, we also used the Northern Hemisphere Equal-Area Scalable Earth Grid 2.0 version 4 weekly snow cover product (NSIDC) to calculate long-term snow disappearance timing and snow cover onset trends since 1980 (Brodzik and Armstrong, 2013; Estilow et al., 2015). The NSIDC product is based on the National Oceanic and Atmospheric Administration (NOAA) climate data record (Robinson et al., 2012). It uses visual interpretation of snow cover detected from a range of sensors (i.e., Advanced Very High Resolution Radiometer – AVHRR, Geostationary Operational Environmental Satellite – GOES, and more recently MODIS; Helfrich et al., 2007) interpolated to the Equal-Area Scalable Earth (EASE) grid with 25 km spatial resolution. The NSIDC product is influenced by image availability and user interpretation of images (Ramsay, 1998; Helfrich et al., 2007). It uses a binary indication of snow or no snow cover. We therefore computed the annual first day with no snow cover for all pixels. Similar to the MODIS product, the snow disappearance timing was determined between 1 February and 31 July. The MODIS and NSIDC snow cover products differ in both their temporal and spatial resolutions, but we found reasonable agreement between snow disappearance timing from both products across the study domain (RMSE = 12.57 Julian day, r=0.76 p<0.01) and individual ecoregions (Fig. S2).

2.3 Fire information

The location and timing of the fire ignitions, as well as their associated burned area, were derived from the Arctic–Boreal Vulnerability Experiment Fire Emission Database (ABoVE-FED) product (Potter et al., 2022). The ABoVE-FED burned area product covers Alaska and Canada (2001–2019) and is derived from thresholding the differenced normalized burn ratio (dNBR) from Landsat imagery at 30 m resolution complemented by MODIS surface reflectance products at 500 m resolution (MOD09GA and MYD09GY v6) when no Landsat data were available. The dNBR thresholding within the ABoVE-FED product was limited to the fire perimeters from the Alaskan Large Fire Database (ALFD, Kasischke et al., 2002) and the Canadian National Fire Database (CNFDB, Stocks et al., 2002), as well as MODIS active fire locations and their surroundings to minimize commission errors from non-fire disturbances (Veraverbeke et al., 2015; Potter et al., 2023).

The retrieval of ignition timing and location was adapted from Scholten et al. (2021b). This algorithm uses the spatiotemporal information in the ABoVE-FED burned area product to delineate individual fire perimeters and a minimum search radius to detect the location of each unique ignition spatially and temporally. Since burned area pixels in boreal regions can be discontinuous due to varying fire severity and possibly omitted pixels, we applied different buffers (1 and 2 km) to group the fire pixels into fire perimeters. Several combinations of the fire perimeter buffers (1 and 2 km), search radii (5, 7.5, 10, and 15 km), and minimum fire sizes (i.e., exclusion of fires from 1 or 2 individual burned pixels) were examined to minimize the commission and omission errors. We tested these three fire size thresholds, as single- or double-pixel burned area could be small anthropogenic fires or commission errors. We compared the results to the ignitions present in the Alaskan Fire Emission Database (AKFED) version 2 (Scholten et al., 2021a) (Table S3). We used ignition locations and timing retrieved inside 2 km buffered fire perimeters with a 7.5 km search radius for fires larger than 50 ha (1 and 2 pixel fires removed) as this was in good agreement with the AKFED-derived ignitions (Table S3). This led to an exclusion of 15 % ignition locations compared to an inclusion of all fire sizes. In Alaska, Yukon, and the Canadian Northwest Territories, we found approximately 6 % more ignitions in ABoVE-FED compared to AKFED and 76 % overlap between the two ignition datasets.

For this study, we also removed ignition locations that were not covered by snow between 2001 and 2019 and ignitions that were erroneously detected before snow disappearance (approximately 11 % of the observations). For the whole study domain and period, we analyzed a total of 17 957 ignitions (Fig. 1b). When possible, we assigned the ignition cause as lightning or anthropogenic from the ignition cause attribute of the ALFD and CNFDB when ignitions fell within the fire perimeter from the same year. By doing so, 4 % of the ignitions were attributed to anthropogenic causes, 38 % were attributed to lightning causes, and the cause of the remaining 58 % was unknown. The daily timing and exact location of fire ignitions were derived from the ABoVE-FED data between 2001 and 2019, but we extended the number of ignitions within ecoregions back to 1980 using fire perimeter data from the ALFD and CNFDB. The start year 1980 was chosen as it corresponds to major optimization of lightning detection systems for Canada that minimized erroneous attribution of causes to fires (Stocks et al., 2002).

We established a relationship between the number of ignitions from ABoVE-FED and the number of fire perimeters from the ALFD and CNFDB for the overlapping period between 2001 and 2019 per ecoregion (Fig. S3). The linear regression between the number of ignitions and fire perimeters was forced through the origin as no fire perimeter can occur without an ignition and vice versa (Fig. S3). The minimum mapping unit (MMU) was 200 ha in CNFDB before 1997 (Stocks et al., 2002) and 405 ha in ALFD before 1988 (French et al., 2015). To minimize uncertainties because of recent changes in the mapping accuracy, we removed fires smaller than 200 ha from the CNFDB and fires smaller than 405 ha from the ALFD, similar to Scholten et al.  (2021b) and Veraverbeke et al. (2017). Similarly, ABoVE-FED fires smaller than MMUs were excluded when developing these relationships. We used the established statistical relationship between ignitions and fire perimeters in each ecoregion from 2001 to 2019 to estimate the annual numbers of ignitions between 1980 and 2000.

2.4 Influence of snow disappearance timing on ignition timing and fire size

For each ignition location, we retrieved the snow disappearance timing by averaging the MODIS-derived day of snow disappearance timing over each ignition location, including its spatial uncertainty derived from the ignition algorithm. Snow disappearance timing may be an important modulator of fire ignitions in the early fire season, whereas seasonal soil moisture dynamics may more importantly influence fire behavior later in the fire season (Flannigan et al., 2016; Gergel et al., 2017). To evaluate the relationship between snow disappearance and ignition timing between 2001 and 2019, we focused on ignitions that occurred early in the fire season. To define early fire season ignitions, we first evaluated the correlation between the annual snow disappearance timing and ignition timing for all ignitions per ecoregion. We then re-evaluated these relationships by only including a fraction of the ignitions. This fraction was derived from taking a percentile of the ignition timing distribution between the 1st and 99th percentile. We generally found significant positive correlations between snow disappearance timing and ignition timing for all percentiles with a general decline in correlation strength with inclusion of ignitions later in the fire season (Fig. S4). Thus, we set the ignition timing threshold to the annual 20th percentile of the ignition timing distribution to account for potential interannual differences in weather and snow disappearance timing interfering with the ignition timing. For this threshold, all ecoregions showed strong significant Pearson r correlation (range: 0.25 to 0.77) between snow disappearance and ignition timing (Fig. S4). By doing so, we retained 3849 ignitions that occurred between the Julian days 58 and 294 across the study domain (Fig. S5).

We also compared all early- versus late-season ignitions to examine the importance of ignition timing for fire size. The burned area caused by an ignition was assigned to the given day of the ignition. In the case of multiple ignition locations detected for one fire perimeter (approximately 4 % of the perimeters), the burned area was assigned to the earliest ignition. We summed up the total burned area between 2001 and 2019 per ignition day. The threshold between early and late ignition timing was again set as the annual 20th percentile day of ignition timing per ecoregion.

2.5 Climatic drivers of snow disappearance and ignition timing

The meteorological drivers of snow disappearance timing and ignition timing were assessed with hourly meteorological data derived from the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis for the climate and weather (ERA5 reanalysis) at 0.25^∘ resolution (Hersbach et al., 2020). ERA5 reanalysis data have been used before in other studies that investigated extreme weather events and fires in the northern high latitudes (Gloege et al., 2022; Parisien et al., 2023). Furthermore, several of the ERA5 variables, such as precipitation, surface temperature, and specific humidity, have been validated with ground observations over the study region (Alves et al., 2020). Fire weather data were collected from the Global ECMWF Fire Forecast ERA5 reanalysis dataset (GEFF-ERA5) of fire danger at 0.25^∘ resolution (Vitolo et al., 2020). We extracted convective potential available energy (CAPE), total precipitation, precipitation type (rain vs. snow), air temperature at 2 m, and dew point temperature at 2 m from the ERA5 reanalysis. From these variables we further derived relative humidity and vapor pressure deficit (Table S4). The fine fuel moisture code (FFMC), duff moisture code (DMC), and drought code (DC) were collected from GEFF-ERA5 and are designed to represent the fuel moisture of the top (1–2 cm organic layer, lag time of $\mathrm{2}/\mathrm{3}$ d), intermediate (5–10 cm sub-organic layer, lag time of 12 d), and deep (15–20 cm deep organic layer, lag time of 52 d) soil layers (Van Wagner, 1987). In regions regularly covered by snow, all fuel load variables are initiated on the third day after the snow has melted, while in regions without snow cover, calculations begin on the third consecutive day with noon air temperatures of <12^∘ (Lawson and Armitage, 2008). Here, we used the fire weather variables as proxies for fuel dryness.

We calculated spatially explicit daily anomalies for all meteorological and fire weather variables by subtracting the climatic daily averages between 1980 and 2019 from the daily observations between 2001 and 2019. We assessed the effect of precipitation, precipitation type (rain vs. snow), air temperature, and relative humidity on snow disappearance timing. Precipitation, air temperature, and relative humidity anomalies were averaged for the 30 d leading up to the day of snow disappearance timing. The number of days with snowfall, rainfall, and no precipitation were summed up for the 30 d leading up to the day of snow disappearance timing. The averages of all weather and fire weather anomalies, excluding precipitation type, between the day of snow disappearance timing and ignition timing were used to assess their influence on ignition timing.

2.6 Temporal trends in snow-free season and fire season

The temporal trends in the snow-free and fire season lengths were analyzed between 2001 and 2019. The snow-free season length was calculated by subtracting the ecoregion average day of snow disappearance timing from the ecoregion average day of snow disappearance offset for each year from the MODIS product.

We evaluated several scenarios to define the fire season timing. For the fire season start, we assessed scenarios between the day of the first ignition and the 20th percentile of the ignition timing distribution. For the fire season end, we assessed scenarios of the day on which 80 % to 99 % of the annual burned area had occurred. First, we analyzed the percentage of annual burned area that was excluded for different fire season start and ending scenarios (Fig. S6). We performed a sensitivity analysis of the different cut-off values that showed no substantial changes in the relationship between the length of the snow-free period and the fire season length (Fig. S7). After evaluation, we chose the 1st percentile day of ignition as fire season start and the day on which 99 % of the annual burned area had occurred as the fire season end day. We subtracted the first day of ignition timing from the day of the 99th percentile total burned area each year to calculate the fire season length. We also investigated changes in the snow-free season length in relation to fire season length between 2001 and 2019.

2.7 Statistical analysis

All statistical analyses were performed in R statistical software version 4.2 (R Core Team, 2022). We investigated temporal trends between 1980 and 2019 and between 2001 and 2019 in snow disappearance timing and the number of ignitions using simple linear regression. The snow-free season length and fire season length in each ecoregion were analyzed between 2001 and 2019 using simple linear regression. The statistical difference in the average fire size between early and late ignitions was analyzed with a Wilcoxon–Mann–Whitney rank-sum test (Mann and Whitney, 1947). We distinguished between two significance levels of p<0.05 and p<0.1.

To assess the ecoregional drivers of the divergent snow disappearance timing and early-season ignition timing, defined as the annual 20th percentile ignitions, we used a linear mixed-effect model. Prior to testing, ignition locations in close proximity were spatially correlated (Moran's I=0.30). We therefore averaged all ignitions for each ecoregion per year to reduce the spatial autocorrelation. The snow disappearance was modeled as a function of weather, while the ignition timing was modeled as functions of weather and fire weather independently. This was to minimize the multi-collinearity in the generalized linear mixed-effect models. We conducted our linear mixed-effect models with ecoregions (16 levels) as random effects using the “nlme” package (Pinheiro et al., 2022) (Tables S5 and S6) (Eq. 1) to account for additional temporal and spatial autocorrelation. We excluded the year as a random effect as it only explained around 3 % and 7 % of the variation in snow disappearance and ignition timing, respectively. We conducted three linear mixed-effect models for all ecoregions combined (n=299), ecoregions with earlier snow disappearance timing trends (n=186), and later snow disappearance timing trends (n=113) based on the MODIS-derived snow disappearance timing (Table S1):

$$y=\mathbf{X}\mathit{\beta }+{\mathbf{Z}}_u+\mathit{\epsilon },$$

where y is the response variable, and Xβ represents the fixed effects, where X is a matrix of observed values per variable and β represents the regression coefficient for each variable. The Z_u term represents the random effects, where Z is a matrix for observed values per covariate of random effects and u is the random effect of the covariates. The error term ε represents the residuals.

All variables were standardized prior to testing and the analysis was conducted for ignitions between 2001 and 2019. The significance of the fixed effects was tested using likelihood ratio tests of the reduced and full models. We used the Akaike information criterion (AIC) to verify the significance of the models compared to reduced models (Zuur et al., 2009). The best model fit was chosen to be the linear mixed-effect model with different intercepts per random effect (ecoregion) bur similar slopes for every predictor and random effect. For further variable selection for our piecewise structural equation model (pSEM), we evaluated the influence of meteorological variables (Table S7) on the day of snow disappearance timing as well as the additional influence of snow disappearance timing and the fuel codes on ignition timing through a redundancy analysis in the R package “vegan” (Oksanen et al., 2013) (Table S7). The significance of the unique contribution of all drivers included in the two variance partitioning analyses was determined by adjusted R² and p<0.05. The shared variance and the residual variance between drivers were also computed (Table S7).

We expected the interactions between predictor variables and the snow disappearance and ignition timing to constitute a complex network and therefore deployed a pSEM in the package “piecewiseSEM” (Lefcheck, 2016). The pSEM creates a single causal network from our deployed multiple linear mixed-effect models that incorporates a random structure (Shipley, 2009). We included explanatory variables linked to snow disappearance and ignition timing based on analysis of bivariate relationships of meteorological and fire weather data that could influence the timing of snow disappearance and ignition. Bivariate relationships were evaluated by simple linear regressions between snow disappearance timing and the respective predictor variables, as well as ignition timing at its potential explanatory variables (Table S8). The hypothesized network of interactions in our pSEM was modeled for three individual pSEMs to test this hypothesized model of interaction between weather, fire weather, and snow disappearance timing but also to describe the potential effect of divergent snow disappearance timing across the study domain. We modeled a pSEM(1) for all ecoregions, (2) ecoregions with early snow disappearance timing trends in accordance with the MODIS trend analysis (Table S1), and (3) ecoregions with later snow disappearance timing trends in accordance with the MODIS trend analysis (Table S1).

For modeling snow disappearance timing, we hypothesized that, (1) as the total amount of precipitation decreases and the air at the surface becomes drier, increased surface air temperature would accelerate snow disappearance timing. We also hypothesized that (2) snow disappearance timing would occur earlier with increased days of no precipitation (smaller snowpack) or days with rain-on-snow events (more rainfall) compared to snow-on-snow events (more snowfall). We hypothesized that (3) earlier snow disappearance timing would result in earlier ignition timing. For modeling the influence of snow disappearance timing on weather and fire weather variables, we hypothesized that (4) surface relative humidity and precipitation would decrease and limit the evaporative cooling, in turn resulting in higher air temperatures. This would increase atmospheric instability, the CAPE, and the likelihood of earlier ignitions. Lastly, we hypothesized that (5) earlier snow disappearance timing would promote drying of fuels (FFMC, DMC, and DC) more pronouncedly in ecoregions with earlier snow disappearance timing. We allow for links between weather and fire weather variables, since DC, DMC, and FFMC are derived from precipitation, relative humidity, and air temperature, while the calculation of FFMC also ingested wind speed. These interactions are included to comply with statistical requirements of inclusion of missing paths in the pSEM analysis but left out of the figure for simplicity reasons (Fig. S8). We also substituted relative humidity and air temperature for vapor pressure deficit in similar pSEMs as VPD has been shown to substantially influence fire ignition and spread (Fig. S9) (Sedano and Randerson, 2014; Veraverbeke et al., 2017). As the pSEMs can consist of many different linear models, we fitted each component of the pSEM with a linear mixed-effect model. Therefore, the influences of fire weather and weather on ignition timing were modeled separately. We included the influence of snow disappearance timing in the model that contained weather variables predicting ignition timing. We assessed potential additional variable interactions and their conditional independence using Shipley's test of dependence separation (d-sep). The test is founded on the χ² distribution and combines the Fisher's C statistics with 2j degrees of freedom, where j is the number of independent interactions in a basis set (Shipley, 2009) (Eq. 2):

$$C=-\mathrm{2}\sum _{i=\mathrm{1}}^k\ln \left(p_i\right)$$

where k is the number of independence claims, and p_i is the null probability of the independence test associated with the ith independence claim.

The missing paths determined by the d-sep test were included in the hypothesized pSEM to accurately analyze the network of dependent variables in our overall pSEM. The global goodness of fit of our models and the hypothesized model was evaluated by the d-sep. With p values > 0.05, the representative model misses no paths and is in accordance with the hypothesized model (Shipley, 2009). The estimates of paths from predictor variables to response variables for each pSEM were standardized for comparison of effects across multiple responses and their indirect and total effects. The standardization of coefficients was done by the ratio of the standard deviation of the independent and dependent variable of the given variables (Eq. 3):

$${\mathit{\beta }}_{\text{SD}}=\mathit{\beta }\times \left({\displaystyle \frac{{\text{SD}}_x}{{\text{SD}}_y}}\right),$$

where β is the unstandardized coefficient, SD_x is the standard deviation of the independent variable, and SD_y is the standard deviation of the dependent variable. The explained variations of snow disappearance and ignition timing from the different components in the pSEMs were analyzed using the marginal and conditional R². Marginal R² represents the variation explained only by the fixed effects, and conditional R² shows the variation explained by a combination of fixed and random effects.

Figure 2 Trends in snow disappearance timing and ignitions for all ecoregions (A–P) between 1980 and 2019 (a). The slope is given for all ecoregions, and its significance level is indicated by ^∗ (p<0.1) or ${}^{\ast \ast }$ (p<0.05). The magnitude and direction of the long-term trends in daily snow disappearance timing (b) and number of ignitions (c) from 1980 to 2019.

3 Results

3.1 Trends in snow disappearance timing and ignitions

The long-term (1980–2019) and short-term (2001–2019) snow disappearance timing trends over boreal North America showed somewhat similar patterns. Long-term snow disappearance timing trends demonstrated shifts towards earlier snow disappearance timing in 13 out of 16 ecoregions, but this trend was only significant in three ecoregions (p<0.05) (Fig. 2a). These significant trends towards earlier snow disappearance were observed in northwestern boreal North America ecoregions (Fig. 2b A–D), while three southern ecoregions (Boreal Plain, Mixed Wood Shield, and eastern Softwood Shield) showed later snow disappearance timing between 1980 and 2019 (Fig. 2a M, O, P). Between 2001 and 2019, this spatial divergence in the trends of snow disappearance timing developed into a distinct west–east divergence across boreal North America. We observed increasingly earlier snow disappearance in western boreal North America versus later snow disappearance in the eastern ecoregions, with only four ecoregions showing statistically significant changes (p<0.05) between 2001 and 2019 (Figs. 1a and 2, and Table S1). The west–east divergence in snow disappearance timing ranged from advances of up to 11 d decade⁻¹ in the western ecoregions to delays of up to 8 d decade⁻¹ in the eastern part of the study region in the MODIS era (2001–2019).

The long-term NSIDC snow product (1980–2019) showed trends between advances in snow disappearance timing of 3 d decade⁻¹ in the west to delays of 2 d decade⁻¹ in the east (Table S1). On average, snow disappearance advanced 1.6 d decade⁻¹ (standard deviation: 0.7 d decade⁻¹, p<0.05) in the western ecoregions (Fig. 2 A–H, J, L), while snow disappearance occurred 1.8 d decade⁻¹ (standard deviation: 0.9 d decade⁻¹) later in the central and eastern ecoregions (p=0.05) (Fig. 2 I, K, M–P). We observed the most pronounced earlier snow disappearance trends of 2.1 d (standard deviation: 0.5 d) earlier snow disappearance decade⁻¹ (p<0.05) in northwestern ecoregions (Fig. 2 A–E), while the most pronounced later snow disappearance trends mainly occurred in the southern ecoregions of 1.1 d decade⁻¹ (standard deviation: 0.8 d decade⁻¹, p=0.05) (Fig. 2 M–P). The spatially diverging trends in snow disappearance timing are associated with similar trends in early spring (February–April) air temperature between 1980 and 2019 (Fig. S10a). The northernmost ecoregions showed the largest increase in early spring air temperature, while the southern ecoregions experienced variation in ignition timing (fire decreasing early spring air temperature over the last 4 decades). Superimposed on this north–south gradient, we also found that the west of the study domain experienced pronounced early spring warming, while the east of the study domain experienced early spring cooling (Fig. S10). The distinct spatial divergence in short-term snow disappearance timing trends also follows a more pronounced short-term early spring air temperature dipole. Early spring air temperatures increased up to 3.5 ^∘C decade⁻¹ over western ecoregions with earlier snow disappearance timing trends and decreased up to 2.1 ^∘C decade⁻¹ over southeastern ecoregions, showing delayed snow disappearance timing (2001–2019) (Fig. S10b).

In accordance with the spatial diverging trends in snow disappearance timing and early spring air temperatures, the trends in the number of ignitions also showed a west–east divergence. The northwestern ecoregions that displayed a pronounced earlier snow disappearance also exhibited an increase in the total number of ignitions of $\mathrm{0.9}\times {\mathrm{10}}^{-\mathrm{6}}$ km⁻² decade⁻¹ (standard deviation: $\mathrm{0.8}\times {\mathrm{10}}^{-\mathrm{6}}$ km⁻² decade⁻¹, p<0.05) (Fig. 2 A–G) between 1980 and 2019. The southwestern ecoregions of the Cold Deserts, Marine West Coast Forest, and Western Cordillera demonstrated the strongest increasing trends in ignitions ($\mathrm{6.4}\times {\mathrm{10}}^{-\mathrm{6}}$, standard deviation: $\mathrm{4.4}\times {\mathrm{10}}^{-\mathrm{6}}$ ignitions km⁻² decade⁻¹, p<0.05) (Fig. 2 H, J, L), while the central and eastern ecoregions showed an overall decrease of $\mathrm{0.2}\times {\mathrm{10}}^{-\mathrm{6}}$ (standard deviation: $\mathrm{0.3}\times {\mathrm{10}}^{-\mathrm{6}})$ ignitions km⁻² decade⁻¹ (p=0.51) (Fig. 2 I, K, M–P). However, there was no spatially divergent trend in the temporal changes in ignition timing between 2001 and 2019. In 12 out of 16 ecoregions, there was a shift towards earlier ignitions when we included all ignitions, with 7 ecoregions showing significantly earlier ignitions (p<0.1). The trends towards earlier ignition ranged between 0.4 and 25 d decade⁻¹ (Table S1). Of the four ecoregions that showed later ignition timing trends, three were located in the southwest of the study domain (Boreal Plains, Cold Deserts, and Western Cordillera), while the western Taiga Shield was the only northern ecoregion that showed a later ignition timing (Fig. S11 and Table S1).

Figure 3 Relationship between snow disappearance and ignition timing for all ignitions of the annual 20th percentile of the ignition timing distribution per ecoregion and their density plots (A–P). The number of ignitions (n), the slope (m), and the significance level (p) are indicated for each ecoregion.

3.2 Relationships between snow disappearance and ignition timing

In all ecoregions, we found significant positive relationships between snow disappearance and ignition timing in the early fire season (20th percentile of the ignition timing distribution) between 2001 and 2019 (p<0.1) (Fig. 3). The strength of the relationships was similar across boreal North America and the advance in ignition timing ranged between 0.22 and 1.43 d per day of earlier snow disappearance (Fig. 3). Ignitions occurred later and in a narrower temporal window in the northern ecoregions (Fig. 3 A–I, K) and eastern Softwood Shield (Fig. 3 P) compared to the other southern ecoregions. Southern ecoregions also showed a more variable ignition timing at the beginning of the fire season (Fig. 3 J, L–P). Furthermore, the southwestern ecoregions of our study domain showed a bimodal ignition timing distribution, which could point to differences in ignition cause. Anthropogenic ignitions dominate earlier in the fire season, while lightning ignitions are more prevalent around the summer solstice (Fig. 3 J, L). Nonetheless, when we separated the anthropogenic and lightning ignitions from ignitions with unknown cause, we still observed positive relationships between snow disappearance and ignition timing for all causes (Table S9).

3.3 Trends in snow-free and fire season lengths

The temporal changes in the snow-free season length and the fire season length also showed a distinct west–east divergence. Corresponding to the overall trends in the snow disappearance timing, we found that the northwestern ecoregions that show increasingly earlier snow disappearance also experience a prolonged snow-free season of 7.1 d decade⁻¹ (standard deviation: 4.2 d decade⁻¹, p<0.1) (Figs. 2a A–H, J, L and 4a A–H, J, L) between 2001 and 2019. The southeastern ecoregions where snow disappearance was occurring later in spring also exhibited a shortening of the snow-free season of 7.3 d decade⁻¹ (standard deviation: 4.7 d decade⁻¹; Figs. 2a I, K, M–P and 4a I, K, M–P), though not significant (p=0.12), between 2001 and 2019. The positive trend in snow-free season length was significant in 5 of the 16 ecoregions, while only the eastern Taiga Shield showed a significant shortening trend in snow-free season length between 2001 and 2019 (p<0.1) (Table S10). We observed similar spatial divergence in the long-term trends in changes in the snow-free season length between 1980 and 2019 (Fig. S12).

The temporal changes in fire season length showed a west–east gradient in complement to a north–south gradient for our study domain (Fig. 4b). The fire season length between 2001 and 2019 increased substantially from 1.7 up to 25.3 d decade⁻¹ and on average 5.8 d decade⁻¹ (standard deviation: 7.6 d decade⁻¹) for the northern ecoregions except in the Taiga Plain (Fig. 4b A–H, K; Table S10) (p=0.45). The southern ecoregions experienced an average shortening of the fire season length between 2001 and 2019 of 18.2 d decade⁻¹ (standard deviation: 10.5 d decade⁻¹) (Fig. 4b I, J, L, M–O) (p<0.1). The northernmost ecoregions in our study region have experienced the largest prolonging of the fire season over the last 2 decades of 18.0 d decade⁻¹ (standard deviation: 10.1 d decade⁻¹) (Fig. 4b B, C, G) (p<0.1).

Figure 4 Changes in the snow-free season length (d decade⁻¹) for all ecoregions (A–P) between 2001 and 2019 (a) and the changes in the fire season length (d decade⁻¹) per ecoregion (A–P) between 2001 and 2019 (b) (Table S10). Letters correspond to the respective ecoregion names (Fig. 2), and significant relationships are indicated by ^∗ (p<0.1) and ${}^{\ast \ast }$ (p<0.05). The relationship between the annual absolute length of the snow-free season from the MODIS product (days) and annual length of the fire season for all ecoregions (c), as well as the ecoregional trends in snow-free and fire season length (d decade⁻¹) (d).

We found that the snow-free season and fire season lengths between 2001 and 2019 were highly correlated (Fig. 4c). There was a consistent significant positive relationship between the snow-free season and fire season lengths across boreal North America between 2001 and 2019 regardless of the thresholds set for the fire season start and end (Fig. S7). Across the study domain, we observed a lengthening of the fire season of 1.7 d for every day of prolonged snow-free season. The length of both the snow-free season and the fire season was shortest in the northern ecoregions and gradually prolonged for more southern ecoregions (Fig. 4b, c). We also found that the trends in snow-free and fire season length tended to correlate positively with each other with a prolonging of the fire season of 0.9 d decade⁻¹ for every day per decade increase in the snow-free season (p=0.08). There was large variation between ecoregions, and the trends in snow-free season lengths explained 45 % of the variation in the trends in fire season length (Fig. 4d).

3.4 Ignition timing and fire size

Fire ignitions that occurred in the early fire season (20th percentile earliest ignitions) resulted in larger fires compared to fires that were ignited later in the season (80th percentile latest ignitions) in all ecoregions but the Alaska tundra (Fig. 5 B) and the eastern Softwood Shield (Fig. 5 P). This difference was significant in 8 out of the 16 ecoregions at p<0.1 with the early-ignited fires resulting in 77 % larger fires compared to fires ignited later in the season across the study domain (Fig. 5). On an ecoregional level, the early-ignited fires grew between 30 % and 600 % larger than late-season fires. The largest difference in fire size between early-ignited and late-ignited fires was observed in the southern ecoregions (Table S10). Also, in these ecoregions, early-season fires accounted for more than half of the total burned area (Fig. 5 J, L, O and Table S10), whereas in the northern ecoregions early-season fires accounted for approximately one-third of the total burned area. Across our study domain, the 20th percentile earliest-ignited fires accounted for an average of 40.6 % (standard deviation: 14.2 %) of the total annual burned area (Table S10). Nonetheless, the largest early-ignited fires on average were observed in the forested ecoregions of the Alaska boreal interior (Fig. 5 A), Taiga Plain (Fig. 5 E), and western Taiga Shield (Fig. 5 G) (23 218 ha, standard deviation: 7557 ha) compared to the other ecoregions (9922 ha, standard deviation: 5192 ha).

Figure 5 Fire size as a function of ignition timing for all ecoregions (A–P). The 20th percentile day of ignition was set as the threshold to discriminate between early- (colored) and late-season fires (gray). The dashed colored lines indicate the mean ignition timing and fire size for all early-season ignitions, while the dashed gray lines indicate the mean ignition timing and fire size for all late-season ignitions. Significant difference between fire sizes of early- and late-ignited fires are indicated by ^∗ (p<0.1) and ${}^{\ast \ast }$ (p<0.05). Note the logarithmic scale for fire size.

3.5 Influence of snow disappearance timing on ignition timing

The pSEM for all ecoregions matched our hypothesized pSEM reasonably well (Fisher's C₈₀=82.24, p=0.41; Fig. 6) and explained 38 % of the variation in the snow disappearance timing (marginal R² (M–R²)=0.38, conditional R² (C–R²)=0.50) and 48 % of the variation in ignition timing (fire weather: M–R²=0.14, C–R²=0.14, and weather: M–R²=0.34, C–R²=0.36) (Fig. 6). The model fits for ecoregions with earlier snow disappearance timing trend (Fisher's C₈₆=96.31, p=0.21) and later snow disappearance timing trends (Fisher's C₁₁₂=107.14, p=0.61) showed similar patterns as the pSEM fit for all ecoregions (Fig. S8). The variances explained in the snow disappearance timing and ignition timing were generally better when splitting ecoregions between those with earlier and later snow disappearance trends. The pSEM for earlier snow disappearance trends explained 32 % of the variation (M–R²=0.32, C–R²=0.32), while 54 % of the variation in ignition timing was explained by the model (fire weather: M–R²=0.15, C–R²=0.15, weather: M–R²=0.39, C–R²=0.39). The pSEM for ecoregions with later snow disappearance trends explained 53 % of the variation in the snow disappearance timing (M–R²=0.53, C–R² = 0.53) and 53 % of the variation in ignition timing (fire weather: M–R²=0.18, C–R²=0.18, weather: M–R²=0.35, C–R²=0.37) (Fig. S8).

These results show that snow disappearance timing was driven by air temperature, without significant influence of precipitation type or amount and humidity. The earlier snow disappearance timing was correlated with high anomalies in air temperature, while the air temperature was generally lower than the climatological averages with later snow disappearance timing (Tables S11–S13). The pSEM results also show that earlier snow disappearance timing promoted fuel drying across ecoregions (Fig. 6 and Table S11). This was also evident from our alternative model, which used vapor pressure deficit instead of relative humidity and air temperature (Fig. S9 and Tables S14–S16).

Snow disappearance timing itself had the strongest individual influence on ignition timing across all ecoregions and models, also after accounting for weather and fire weather. The cascading effect of accelerated drying of organic soils from earlier snow disappearance timing carried over to the timing of ignition. For all models, the DMC had the strongest influence on the ignition timing, while the FFMC significantly affected ignition timing across all ecoregions and over the ecoregions exhibiting earlier snow disappearance timing (Figs. 6 and S8a). For ecoregions with later snow disappearance trends, only the slow responding fuel moisture codes (DMC and DC) significantly influenced the timing of ignition. For ecoregions with earlier snow disappearance timing, DC influenced the ignition timing positively, meaning that earlier ignitions generally occurred when DC was still low. The fuel moisture codes together influenced the ignition timing more strongly compared to snow disappearance timing and weather variables across models (Tables S11–S13).

Early snow disappearance may also affect larger-scale atmospheric dynamics. We found that earlier snow disappearance timing contributed to the destabilization of the atmosphere through increased CAPE across ecoregions (Fig. 6), in particular for ecoregions with earlier snow disappearance timing (Fig. S8a). Early snow disappearance was associated with higher air temperatures and lower humidity in the overall model. These favorable weather conditions led to earlier ignition in the overall model and the model for ecoregions with earlier snow disappearance timing. Early ignitions were associated with lower relative humidity and higher air temperatures driven by the earlier snow disappearance timing (Figs. 6 and S8a). Snow disappearance timing itself had the strongest individual influence on ignition timing across all ecoregions and models. Similar results were obtained when air temperature and relative humidity were substituted by vapor pressure deficit (Fig. S9 and Tables S14–S16).

4 Discussion

4.1 Diverging spatial trends in snow disappearance timing and ignitions

We found the co-occurrence of a pronounced continental dipole in decadal trends of snow disappearance timing and number of fire ignitions across arctic–boreal North America. We observed increasingly earlier snow disappearance and an increase in the number of fire ignitions in northwestern boreal North America between 1980 and 2019. In contrast, snow disappearance timing has simultaneously been occurring later and the number of fire ignitions decreased in the last decades in the southeastern part of our study domain. The divergent trend in the number of ignitions is in accordance with a previous study on changes in the number of fires and amount of burned area across Canada from 1959 to 2015 (Hanes et al., 2019). The changes in snow disappearance timing that we found in our study are corroborated by earlier work demonstrating both increasing and decreasing trends in snow cover over southeastern and northwestern boreal North America, respectively (Chen et al., 2016). Furthermore, Bormann et al. (2018) found an earlier onset of spring snow disappearance in northwestern boreal North America in contrast to later snow disappearance or no changes in snow disappearance timing over southeastern boreal North America.

We also found that the west–east diverging trend in snow disappearance timing has become more pronounced in the last 2 decades compared to the longer-term trend since 1980. These observations follow the divergent trend of less pronounced changes in long-term early spring air temperature (1980–2019) and distinct dipoles in early spring air temperature over boreal North America between 2001 and 2019 (Fig. S10). Similar to Cohen et al. (2014), we found small changes in air temperature between 1980 and 2019 in the northern and southern parts of our study domain (Fig. S10a). The last 2 decades of enhanced west–east divergence in snow disappearance timing followed the development of a pronounced west–east dipole of trends in early spring air temperature as observed in our linear mixed-effect models (Table S5) and also observed in two consecutive recent winters between 2013 and 2015 (Singh et al., 2016). As higher early spring air temperatures promote earlier snow disappearance, the snow–albedo feedback will in turn result in higher air temperatures (Déry and Brown, 2007). In this way, the presence of a dipole of changes in early spring air temperature and snow disappearance timing over boreal North America might indicate that both processes reinforce each other on subcontinental scales.

Besides regional changes in early spring air temperature, large-scale atmospheric dynamics may have also influenced snow disappearance timing and the number of ignitions as observed in our study (Cohen et al., 2014; Zhao et al., 2022). Changes in sea ice and snow cover (Zou et al., 2021) may have a large impact on the location of the polar jet stream and tropospheric ridge persistency, causing air temperature extremes (Francis and Vavrus, 2012; Kim et al., 2014; Horton et al., 2016). In recent decades, these persistent tropospheric ridge patterns have been located over the northwestern part of our study domain, which traps and slows the progression of Rossby waves eastwards (Francis and Vavrus, 2012; Jain and Flannigan, 2021), resulting in downstream troughing over the east (Singh et al., 2016). This tropospheric ridge leads to a blocked anticyclone in the west, causing higher air temperatures and increased burned area, and an associated cyclone over eastern North America with lower air temperatures and less burned area (Skinner et al., 1999; Cohen et al., 2014; Sharma et al., 2022). Further, the stratospheric vortex, which is made up of westerly winds formed in the stratosphere during wintertime, may have weakened, and consequently sudden stratospheric warming (SSW), which is rapid heating in the stratosphere over the North Pole, has caused winter cold spells over eastern Canada over the last 4 decades (Kretschmer et al., 2018b). The effect of winter cold spells may carry over into spring, delaying the snow disappearance timing and thus the fire season.

The presence of a dipole in snow disappearance timing and ignition trends in our study is likely related to (1) changes in the stratospheric vortex and SSW that send winter cold spells over the eastern part of the study domain (Kretschmer et al., 2018a) and as a consequence annual mean air temperature anomalies divergence from increasing in the west to decreasing in the east of boreal North America in the last decades (Cohen et al., 2012; Coumou et al., 2018) as well as (2) changes in the location of the summer jet as a consequence of longer persistence of positive geopotential anomalies over the western part of our study domain (Jain and Flannigan, 2021; Zou et al., 2021). Although our results do not provide clear indications, the persistent ridge formation could possibly also be a result of the divergent snow disappearance trend caused by the SSW. Both the soil moisture and albedo feedback between snow disappearance timing and air temperature may have further strengthened the diverging trends. These processes may have also influenced the fire extremes across Canada in 2023. Our results suggest that these atmospheric processes and soil moisture feedbacks may have also led to the enhanced fuel dryness in western ecoregions that has driven the large increases in the number of ignitions compared to the other ecoregions (Abatzoglou and Williams, 2016; Holden et al., 2018).

4.2 Influence of snow disappearance timing on ignition timing and fire size

By focusing on the start of the fire season, we were able to disentangle the effect of snow disappearance timing on ignition timing. Previous studies found no significant effects of snow disappearance timing on annual burned area, with snow disappearance timing being regarded as a minor driver of annual burned area compared to meteorological variables (Jolly et al., 2015; Kitzberger et al., 2017). Nonetheless, snow disappearance timing has been shown to play a crucial role in altering fuel dryness and the frequency of large fires over a temperate forest in the western United States (McCammon, 1976; Westerling et al., 2006). Our results show that snow disappearance timing has a strong influence on early ignition timing and fire size in all ecoregions of boreal North America. This relationship diminished when snow disappearance timing was compared to progressively later ignitions (Fig. S4). This may be due to the importance of the spring window, the period between snow disappearance timing and leaf flush, on early-season fires (Parisien et al., 2023). During the spring window deciduous and mixed forests are very conductive to fire and ecoregions experiencing the longest spring window corresponded to where we also found the highest early fire ignition density (Fig. 5 J, L–M) (Parisien et al., 2023). Also, the longest spring window was found in the interior west of Canada (Parisien et al., 2023), which coincides with the ecoregions with most fire ignitions observed in our study (Fig. 1). Late-season ignitions in July, August, and September may be more influenced by long-term drought and synoptic weather conditions than by snow disappearance timing (Jain et al., 2017; Holden et al., 2018).

We found that throughout boreal North America, fires caused by early-season ignitions following earlier snow disappearance also on average grew larger than fires ignited in the late fire season. This was in accordance with earlier findings limited to Alaska, USA, and the Canadian Northwest Territories (Veraverbeke et al., 2017). Because of the early snow disappearance and the earlier ignition timing, early-season fires have longer temporal windows with potential for favorable warm and dry weather conditions conductive to fire spread (Sedano and Randerson, 2014). Indeed, the 20 % earliest ignitions resulted in approximately 40 % of the total burned area across the study domain between 2001 and 2019. In the future, the contribution of early fires to burned area might increase with warmer and drier weather conditions, leading to earlier snow disappearance and thus an increased likelihood for earlier and larger fires over boreal North America (Flannigan et al., 2005, 2013).

Figure 6 Piecewise structural equation model (pSEM) of the hypothesized snow disappearance and ignition timing model (a) and its fit for all ecoregions (b). The influences of fire weather and weather on ignition timing were modeled separately. Gray arrows represent positive effects and black arrows indicate negative effects. The single-headed arrows show significant direction of causal relationships, while double-headed arrows represent significant non-causal relationships (p<0.1). All arrows are scaled to their respective effect size (Table S11). Marginal R² (M–R²) indicates the variation solely explained by the fixed effects, and conditional R² (C–R²) represents the variation explained by both the fixed and random effects.

4.3 Changes in the snow-free and fire season lengths

We found a north–south gradient in the changes in the actual fire season length ranging from a prolonging of 30 d decade⁻¹ in northern ecoregions to a shortening of 25 d decade⁻¹ in southern ecoregions. Previous studies have mainly found the prolonging of the potential fire season to be between 3 and 30 d decade⁻¹ over boreal North America (Wotton and Flannigan, 1993; Flannigan et al., 2013; Jolly et al., 2015; Jain et al., 2017). These estimates of the prolonging of the potential fire season were based on changes in fire weather (Flannigan et al., 2013; Jain et al., 2017). Other studies that examined governmental fire perimeter data also found a prolonging of the fire season length limited to the western North America (Westerling et al., 2006; Albert-Green et al., 2013; Hanes et al., 2019). In our study, we used daily fire spread data derived from satellite observations to determine the fire season start and end dates (Skakun et al., 2021; Potter et al., 2023). This approach, however, relies on MODIS active fire data and is therefore limited to the MODIS era in the 2000s. Longer-term accurate temporal and spatial data on ignition timing and end of burning are needed to assess the changes in the actual fire season on a climatic timescale and a continental scale. Our results suggest that a change in the duration of the snow-free season is almost one-to-one related to a change in the duration of the fire season across boreal and Arctic North America. However, the effect may be of more importance for the fire season start than the end of the fire season as this is often marked by the first rainfall in autumn in adequate amounts to extinguish fires and rewet dried-out fuels, preventing new ignitions. Climate-change-induced changes in the amount and timing of autumn rainfall will likely affect the timing of the fire season end (Holden et al., 2018; Goss et al., 2020), although recent studies also showed that some fires overwinter and re-emerge the following spring (McCarty et al., 2020; Scholten et al., 2021b; Xu et al., 2022), challenging the concept of a demarcated fire season. In a warmer North American arctic–boreal region, the snow-free season will likely be prolonged with a consequent lengthening of the fire season, both starting earlier in spring and being prolonged later into autumn (Flannigan et al., 2013).

4.4 Cascading effects of snow disappearance timing on weather and ignition timing

In the three piecewise structural equation models (pSEMs), anomalies in snow disappearance timing were only attributed to anomalies in air temperature (hypothesis 1). Our models did not confirm the importance of the amount or the type of precipitation for snow disappearance timing observed in previous research (hypothesis 2) (Barnett et al., 2005; McCabe et al., 2007). However, air temperature also affected precipitation types in our models, which, although statistically insignificant, showed divergent influences on snow disappearance timing between ecoregions with earlier and later snow disappearance trends (Tables S12 and S13). This suggests that the air temperature dipole observed in the last 2 decades (Fig. S10) may influence precipitation, including snowpack volume and persistency (Brown and Mote, 2009) and therefore likely also snow disappearance timing (Barnett et al., 2005). Nonetheless, snow disappearance timing itself largely influenced ignition timing regardless of ignition source and ecoregion, and we additionally found snow disappearance timing to be an important early indicator for early-season fires in the North American arctic–boreal region (hypothesis 3). Our model also suggests a cascading effect of snow disappearance timing on meteorological conditions that carried over into the influence on ignition timing. Relationships between warm and dry conditions and ignitions and fire spread have been established before (Sedano and Randerson, 2014; Veraverbeke et al., 2017). This was only apparent for the overall model and the model including ecoregions with earlier snow disappearance timing (hypothesis 4). This suggests that land–atmosphere dynamics are altered by changes in snow disappearance timing as it influences the soil moisture content, which is proportional to evapotranspiration, changing the land energy balance (Seneviratne et al., 2010). Also, this in combination with anomalously high springtime air temperatures promoted greening of vegetation and desiccation of soils in other boreal regions, changing the impact of the warming in the atmospheric (Gloege et al., 2022). These land–atmosphere dynamics may have been potential pathways for extreme fire seasons in Siberia (Scholten et al., 2022), and our pSEM results suggest that similar dynamics may be in place over ecoregions with earlier snow disappearance timing. The ecoregions with later snow disappearance timing, which showed less of a carry-over effect of snow disappearance timing on weather and fuel moisture to ignition, also corresponded to the more densely populated regions with larger fractions of deciduous trees (Pavlic et al., 2007; Olthof et al., 2015). The lack of a carry-over effect may be due to other drivers such as elevated potential for anthropogenic ignitions and the spring window during which deciduous forests are most flammable (Wotton et al., 2010; Parisien et al., 2023).

The results of our study also point to cascading effects of changes in snow disappearance timing on dry fuel availability that carried over into the ignition timing across all models. The fine fuel moisture and duff moisture codes showed significant influences on ignition timing, while the drought code did not, possibly due to its slower response to changes in weather (hypothesis 5). This is in agreement with previous studies that indicate that the ignition of fires in boreal North America strongly depends on the immediate dryness of the fine fuels (Abatzoglou and Williams, 2016; Hessilt et al., 2022). The effects of earlier snow disappearance timing on enhanced desiccation of fuels observed in three forests sites (McCammon, 1976) may be broadly applicable across boreal North America. As observed in our study, dry fuels can directly promote ignition timing as they are readily ignitable (Hessilt et al., 2022), but this may also be indirect through the influence on aboveground biomass senescence and ecosystem production (Liu et al., 2020). Lastly, the changes in fuel moisture as a consequence of snow disappearance can also play a role in the soil moisture–climate interaction, which is not accounted for in this model.

4.5 Limitations

We used a conservative threshold of 14 consecutive days of snow-free pixels (NDSI ≤ 15 %) to calculate the snow disappearance timing. This could potentially cause the timing of snow disappearance to occur later than observed. A comparison of snow disappearance timing retrieval with a threshold of 7 consecutive days of snow-free pixels indicated that the retrievals resulted in similar temporal patterns in snow disappearance timing regardless of threshold choice (Fig. S1, Table S2), with the 14 d threshold generally resulting in later snow disappearance timing. The largest discrepancies between the retrievals of snow disappearance timing with different temporal thresholds were found in the southern ecoregions (Table S2). This indicates that the threshold of 14 consecutive days with snow-free pixels may be more robust to determine snow disappearance timing because sudden changes in weather can manifest in snow offset and onset, especially in southern ecoregions.

The long-term and short-term trends of snow disappearance timing and the number of ignitions were not consistently statistically significant for all ecoregions. The uncertainty related to the retrieval of both snow disappearance timing and fire perimeters from the pre-MODIS era (before 2000) resulted in large variation in both variables. More robust findings could potentially be obtained with longer time series. Also, continued monitoring of snow disappearance and ignition timing is needed to track the relationship between these two variables as their relationship may become more pronounced with further climate change. Similarly, longer and more consistent time series could increase the robustness of the analysis of snow-free vs. fire season length. While we observed a significant relationship between these variables across ecoregions (Fig. 4c), this was not evident for most ecoregions (snow-free periods: p<0.1 in six ecoregions and fire season length: p<0.1 for two ecoregions). We only used the shorter time period between 2001 and 2019 data to establish these changes as these represent higher-quality data during the MODIS era.

Lastly, our pSEM analysis gives a simplified overview of relationships between snow disappearance timing, land–atmosphere dynamics, and fire ignitions. However, we acknowledge that these interactions are highly coupled. The complexity is beyond our model and may involve variables that we did not include. The influence of snow disappearance timing on atmospheric variables through surface albedo change and altered soil moisture may be difficult to decouple from the atmospheric variables and their persistent seasonal patterns in snow disappearance timing itself. Our models provide further support of the importance of land–atmosphere dynamics in relation to fire, yet our analysis did not provide robust relationships explaining the mechanistic interactions. We therefore call for a better understanding of the role of snow disappearance timing in land–atmospheric dynamics affecting boreal fires, specifically large-scale influence of continuous snow disappearance on soil and fuel properties, e.g., soil and fuel moisture, as well as atmospheric conditions, e.g., vapor pressure deficit, and vice versa. Understanding these interactions and feedbacks could further advance our comprehension of how climate change is affecting changing boreal fire regimes.

5 Conclusions

We found a pronounced west–east divergence of recent changes in snow disappearance timing and the number of fire ignitions across boreal North America. Our results point to a clear trend of earlier spring snow disappearance in the northwestern ecoregions, while the southern and eastern ecoregions showed an increasingly later snow disappearance timing over the last decades. Similarly, the total number of fire ignitions increased in the northern and western ecoregions, while the southeastern ecoregions experienced little to no changes in the number of early fire ignitions. We conclude that climate warming resulted in increasingly earlier snow disappearance in northwestern boreal North America, which in turn led to earlier fire ignitions that tended to grow into larger fires.

The temporal trends in snow disappearance and ignitions across boreal North America followed the same spatial pattern of temporal trends in early spring air temperature over the last 4 decades. Snow disappearance and ignition timing were positively correlated across all ecoregions, and earlier snow disappearance was also the main driver of earlier fire ignitions across all ecoregions. Further, we found a cascading effect of elevated air temperature and earlier snow disappearance that carried over into earlier drying of fuels, which resulted in earlier ignitions across the study domain. This cascade was more pronounced over ecoregions with increasingly earlier snow disappearance timing than over those with increasingly later snow disappearance timing. Our work points to the important impact that snow cover and snow disappearance timing have on fire ignitions and fire size across boreal North America, as well as the influence of changes in snow disappearance timing on changes in fire regimes. In a warming North American boreal forest, earlier snow disappearance will likely result in increasingly earlier and larger fires.