Observed links between heatwaves and wildfires across Northern high latitudes

This study focuses on the terrestrial ASA (60°–80°N, 180°W–180°E). We divide the study area into quadrants: 180°–90°W (Q1), 90°–0°W (Q3), 0°–90°E (Q4), 90°–180°E (Q2). The choice of quadrants follows from previous research by Dobricic et al (2020). Due to the short warm period in the ASA, we limit our HW analysis to June–August (JJA) and BA analysis to June–September (JJAS). We restrict our analysis to 2000–2022 due to data availability. We regrid each dataset to 1$^\circ\times$1° resolution to identify large-scale heat and fire patterns.

We define burnable area within our study region using the Boreal-Arctic Wetland and Lake Dataset (BAWLD). It provides estimates of fractional land cover of 19 land cover classes at 0.5° resolution (Olefeldt et al 2021). The study region of BAWLD is categorised into 15 wetscapes via distinct configurations of land cover classes. We mask out wetscapes classed as Rivers, Large Lakes, Glaciers, Alpine and Tundra Barrens and Lake-rich Shield, as they comprise of little or no vegetation cover. We assume the other ten wetscapes have sufficient vegetation that could burn—we call these regions burnable.

We use the LST/Emissivity Daily (MOD11C1) Version 6.1 product (Wan et al 2021) from moderate resolution imaging spectroradiometer (MODIS) aboard Terra from 2000 to 2022. MOD11C1 provides gridded, daily LST at 0.05° resolution. The average overpass time of Terra is around noon. As the LST algorithm only works for clear-sky pixels, cloud cover leads to missing values. This means that MOD11C1 has an intrinsic warm bias in its climatology; it under-represents cloudy and rainy days, which tend to be cooler. However, MOD11C1 does not under-sample clear, warm days. Heavy smoke might also hinder observations. Comparison of LST and surface air temperature from ERA5 shows strong correlation between the data, while confirming an LST warm bias (see supplementary materials). While several studies attempt to create global, spatio-temporally continuous LST datasets (Ghafarian Malamiri et al 2018, Yu et al 2022, Zhang et al 2022), we do not interpolate missing values in order to maintain purely observational data.

We retrieve two fire products derived from MODIS:

• Monthly BA from MODIS (FireCCI51), produced as part of the Fire Disturbance Climate Change Initiative project of the European Space Agency (Chuvieco et al 2018). FireCCI51 is available for 2001–2020 at 250 m resolution, we aggregate it to 1$^\circ\times1^\circ$ to focus on large fires.
• Archived active fire data (MCD14ML) distributed by Fire Information for Resource Management Service, spans 2001–2022 (NASA FIRMS 2023). MCD14ML is a subset of the data produced by the MODIS Fire and Thermal Anomalies algorithm (Giglio et al 2016), aboard both Terra and Aqua. MCD14ML records $1\times1\,\mathrm{km}$ pixels containing one or more active fires or thermal anomalies. We aggregate the data to 1$^\circ\times1^\circ$ resolution by summing up the number of active fire pixels within each grid cell.

Neither fire product is available in 2000, so HWs in 2000 are excluded from our HW–wildfire analysis. As a level 2 product, MCD14ML has spatial and temporal sampling issues (Giglio et al 2016) in the presence of smoke. Hence the monthly BA product may be more robust and representative of wildfire extent on seasonal timescales. We use MCD14ML to study the relationship between HWs and wildfires on timescales shorter than a month and for investigating the timing of peak fire activity with respect to peak HW activity.

Our HW definition builds on previous definitions used in climate research (Perkins-Kirkpatrick and Alexander 2013, Cowan et al 2016). On a grid-point level, we define heatwaves as three consecutive days which exceed the day-of-year 90th percentile (90thp) of LST. We use a multi-year, centred 11 day rolling window for calculating 90thp for each grid point over the analysis period. If a missing day caused by sensor outage in 2000 is both preceded and followed by an exceedance, we flag it as an HW pixel. We use the 90thp, instead of the 95th percentile, due to the presence of missing values. We use a rolling window of 11 days to account for the seasonal cycle, provide a reasonable sample size and reduce noise from synoptic scale variability. We use at least three consecutive days of anomalous heat, as more consecutive days yield more intense events and prolonged heat is more likely to impact vegetation. Reasonable perturbations of the percentile and length of consecutive days to identify gridpoint heat, and the length of the rolling window for climatology yield similar results.

To identify space-time coherent, large-scale HWs, we perform spatio-temporal density-based clustering with noise (ST-DBSCAN) on our dataset of grid points recording HW conditions. We apply a three-dimensional (3D) version (Birant and Kut 2007) of DBSCAN (Ester et al 1996), an unsupervised clustering algorithm that labels each point as part of a group or noise based on its location with respect to all other points in the dataset. We discard small clusters below 100 000 km² daily extent (details see supplementary materials) to study the characteristics of large HWs in the Arctic. This results in the final definition of a (large) HW as prolonged heat exceedances clustered in space and time, which cover at least 100 000 km² for at least one day.

The 3D ST-DBSCAN requires three parameters (Birant and Kut 2007, Cakmak et al 2021): a spatial threshold ($\varepsilon_\mathrm{s}$), the minimum number of neighbours (η) and a temporal threshold ($\varepsilon_\mathrm{t}$). $\varepsilon_\mathrm{s}$ and $\varepsilon_\mathrm{t}$ are the maximum spatial and temporal separation between two points to be neighbours and η is the minimum number of points within a neighbourhood for the region to be considered dense.

We evaluate the sensitivity of HW clusters to parameter choices. They are least sensitive to $\varepsilon_\mathrm{s}$, but vary greatly with $\varepsilon_\mathrm{t}$ and η. The clustering performed best when choosing 1000 km, 2 days and 10 neighbours as parameters; see supplementary materials for a detailed description of our parameter tuning method and discussion on clustering algorithm choices. We also note that the clustering algorithm is a tool to identify larger and longer HWs, rather than find the optimal HW characteristics. We compare clusters found in the LST data to those in the daily maximum near-surface air temperature (${T_{\text{max}}}$) from ERA5 reanalysis (table S2). We find good agreement, but more clusters in ${T_{\text{max}}}$ due to the aforementioned warm bias in LST data (see supplementary materials).

We are not aware of any other study that uses ST-DBSCAN for HW analysis. Most HW-related studies define HWs at grid cell level only (Reddy et al 2022). To characterise HW extent, many studies (e.g. Cowan et al 2016) estimate the proportion of grid cells under HW conditions within a region of interest, instead of tracking the shape and extent of the events.

Using only grid points which are a part of an HW cluster without interpolation, we assign attributes to each cluster:

• Duration: the number of days from the first to last day when the daily cluster size exceeds 100 000 km². This filters out small HWs, the beginning and tail of large HWs. However, it lets HWs vary in size throughout their duration.
• Region: the quadrant that the coordinates of the cluster centroid are located in.
• Amplitude: maximum of the daily area-weighted mean LST anomaly of the hottest 100 000 km² within a cluster. For each day, we rank the highest LST anomalies from the 90thp and cut off our weighted average once the cumulative extent of the ranked grid cells exceeds 100 000 km². This method samples over the hottest regions of the same size to avoid a low bias in the intensity of the largest HWs.
• Entropy: scatter of LST anomalies within an HW. Higher entropy means more variation in LST anomalies.
• Anomalous energy: integral of the daily LST anomaly over the spatial extent of the HW, averaged over its duration.
• Maximum extent: HW area on the day of the largest spatial extent.
• Fire density: mean fire count per day and per burnable area. We divide the number of active fires during each day of a HW by its burnable area, and average over its duration. The resulting fire density values are standardised by the spatial mean and standard deviation (σ) in each quadrant separately—0 corresponds to the regional-average fire activity during the HW, and 1 (or above) means that the fire occurrence was 1 σ (or more) above the regional-average.

To compare HWs within and across the quadrants, we use the Pearson correlation coefficient (r) and the multiple correlation coefficient ($r_\mathrm{m}$, for $\unicode{x2A7E}3$ variables).

We aggregate HW days over summer (JJA) as the percentage of burnable grid cells with valid LST observations which are part of a cluster. Similarly, we divide the cumulative sum of BA over the extended summer (JJAS) by its total burnable area. We estimate the Pearson correlation coefficient r for the resulting time series in each quadrant. This gives an indication of the strength of the HW–wildfire relationship. We assess the significance of r by randomly perturbing the time-order of the aggregated BA, with p-value being the proportion of 10 000 permutations which has a larger correlation coefficient than the observed r.

We also analyse the spatial characteristics of the mean JJAS BA in each quadrant. We compare these results to spatial distribution of the fire density of HW clusters during the whole study period, the five strongest HW and fire summers.

To analyse timing between heat and fire, we create composites of HW clusters with the highest amplitude and stack them by the day of their maximum amplitude. We carry out lagged cross-correlation (CC) of the daily HW metric and fire count. We perform a leave-one-out cross-validation to test sensitivity of our composite mean to each member.

The HW attributes and number of clusters we obtain are shown in figure 1. We note no significant difference between the distributions of HW characteristics across the four quadrants (figure 1(a)). We find that HW metrics are significantly correlated (p < 0.05) with each other in all quadrants: max extent with duration $\mathrm{(Q1:}\;r\! = \!0.65\!\,\pm\!\,0.17\mathrm{,\;Q2:\;}r\! = \!0.59\!\,\pm\!\,0.16\mathrm{,\;Q3:\;}r\! = \!0.51\!\,\pm\!\,0.25\mathrm{\;and\;Q4:\;}r\! = \!0.61\!\,\pm\!\,0.17\mathrm{)}$, max extent with amplitude $\mathrm{(Q1:}\;r\! = \!0.63\!\,\pm\!\,0.17\mathrm{,\;Q2:\;}r\! = \!0.70\!\,\pm\!\,0.13\mathrm{,\;Q3:\;}r\! = \!0.76\!\,\pm\!\,0.15\mathrm{\;and\;Q4:\;}r\! = \!0.69\!\,\pm\!\,0.14\mathrm{)}$, and duration with amplitude $\mathrm{(Q1:}\;r\! = \!0.44\!\,\pm\!\,0.23\mathrm{,\;Q2:\;}r\! = \!0.48\!\,\pm\!\,0.19\mathrm{,\;Q3:\;}r\! = \!0.30\!\,\pm\!\,0.31\mathrm{\;and\;Q4:\;}r\! = \!0.61\!\,\pm\!\,0.17\mathrm{)}$. In addition, larger and longer HWs have higher amplitude $\mathrm{(}r_\mathrm{m}\! = \!0.63\mathrm{\;in\;Q1,\;}r_\mathrm{m}\! = \!0.70\mathrm{\;in\;Q2,\;}r_\mathrm{m}\! = \!0.77\mathrm{\;in\;Q3,\;}r_\mathrm{m}\! = \!0.73\mathrm{\;in\;Q4)}$. Figure 1(b) demonstrates these relationships visually for each quadrant.

Zoom In Zoom Out Reset image size

With the active fire archive, we also find significant correlation between fire density and maximum extent in all quadrants $\mathrm{(Q1:}\;r\! = \!0.37\!\,\pm\!\,0.24\mathrm{,\;Q2:\;}r\! = \!0.26\!\,\pm\!\,0.23\mathrm{,\;Q3:\;}r\! = \!0.38\!\,\pm\!\,0.29\mathrm{\;and\;Q4:\;}r\! = \!0.51\!\,\pm\!\,0.20\mathrm{)}$, and with other characteristics in some quadrants $\mathrm{(for\;duration:\;Q1:}\;r\! = \!0.22\!\,\pm\!\,0.23\mathrm{,\;Q2:\;}r\! = \!0.56\!\,\pm\!\,0.17\mathrm{,\;Q3:\;}r\! = \!0.14\!\,\pm\!\,0.33\mathrm{,\;Q4:\;}r\! = \!0.46\!\,\pm\!\,0.21\mathrm{;\;and\;for\;amplitude\;Q1:\;}$ $r\! = \!0.19\!\,\pm\!\,0.24\mathrm{,\;Q3:\;}r\! = \!0.20\!\,\pm\!\,0.32\mathrm{,\;Q4:\;}r\! = \!0.21\!\,\pm\!\,0.25$) despite substantial sampling uncertainty. The strong correlation between fire density and maximum extent is consistent with drying of the vegetation caused by large-scale extreme heat promoting burning (Seneviratne et al 2021), and increased chance of ignition in a larger area.

We also identify a link between HW occurrence and wildfire BA across a quadrant over the whole summer (figure 2). While the relationship between JJA HW and JJAS BA remains significant over the entire study region (r = 0.63), splitting the Arctic into quadrants reveals that the quadrants containing Siberia (r = 0.87) and Alaska (r = 0.46) are the strongest contributors toward the observed strong HW–BA relationship within the Arctic. Q3 and Q4 show low, non-significant correlation between summertime HW activity and BA (r = 0.04 and r = 0.27 respectively). This is consistent with less vegetation cover in NE Canada and Greenland (Q3). The reason for this smaller correlation for Q4 (encompassing North Europe) is less clear, and may be linked to a strong fire management regime in North Europe (Nordic Forest Research, SNS 2021) and possibly also complex topography. As the Arctic warms, HWs will become more likely (Seneviratne et al 2021), but due to high variability this occurs on the background of large climate variability seen in the sectors. For discussion on trends and the impact of clustering on figure 2, see supplementary materials.

Zoom In Zoom Out Reset image size

We further investigate the temporal link between HW and fire activity in figure 3, using MCD14ML. As the quadrant (Q2) exhibiting the strongest HW–wildfire link, we create a composite of the ten hottest HW clusters in Siberia in figure 3 (see supplementary materials for other quadrants). Peak fire activity tends to follow peak HW activity—significant (p < 0.05) CC is found for a lag of 12–13 d. While the location of significant heat/fire lags varies with the cluster number or HW metric chosen, the strong and significant negative correlation before the peak HW activity remains consistent, showing that fire activity does not climax before HW activity.

Zoom In Zoom Out Reset image size

We investigate the spatial relationship between HW and fire activity in figure 4—relating fire activity (BA as grayscale shading and fire density as coloured points) to HW locations. Evaluating the entire study period (figure 4(a)), we find negligible fire activity during HWs in Greenland and North Europe. Events coinciding with anomalously strong fire activity (in space, at any given time) in Q4 are located in Western Siberia. In Q1 and Q2 we see some overlap between where high fire density HWs have been observed and where most of the summertime BA occurs, and a tendency for anomalous fire density (red shaded dots) to occur with larger HW (larger dots) in those regions. This indicates a spatial dependency for fire-prone HWs in Interior Alaska, the Canadian Prairies and Eastern Siberia. While HWs occur across the entire Arctic, there is a limited relationship between HWs and fires in the far northern Arctic. In summary, extreme heat over regions with sufficient vegetation appear most likely to result in high fire activity.

Zoom In Zoom Out Reset image size

Finally, we investigate whether HW activity in the five strongest HW/fire summers is more prone to fires (figures 4(b) and (c)). In comparison to figure 4(a), we find a more pronounced spatial HW–BA relationship during these five summers (figures 4(b) and (c)). Especially for intense fire summers, we find a strong spatial co-incidence between HWs and BA (figure 4(c)). This suggests that extensive HWs contribute greatly to intense fire summers. As fire activity lags behind peak HW activity (figure 3), this supports that intense HW activity may result in extreme BA. However, strong HW activity does not result in extensive burning everywhere (figure 4(b), and discussion).