Annual and Seasonal Patterns of Burned Area Products in Arctic-Boreal North America and Russia for 2001–2020

: Boreal and Arctic regions have warmed up to four times quicker than the rest of the planet since the 1970s. As a result, boreal and tundra ecosystems are experiencing more frequent and higher intensity extreme weather events and disturbances, such as wildfires. Yet limitations in ground and satellite data across the Arctic and boreal regions have challenged efforts to track these disturbances at regional scales. In order to effectively monitor the progression and extent of wildfires in the Arctic-boreal zone, it is essential to determine whether burned area (BA) products are accurate representations of BA. Here, we use 12 different datasets together with MODIS active fire data to determine the total yearly BA and seasonal patterns of fires in Arctic-boreal North America and Russia for the years 2001–2020. We found relatively little variability between the datasets in North America, both in terms of total BA and seasonality, with an average BA of 2.55 ± 1.24 (standard deviation) Mha/year for our analysis period, the majority (ca. 41%) of which occurs in July. In contrast, in Russia, there are large disparities between the products—GFED5 produces over four times more BA than GFED4s in southern Siberia. These disparities occur due to the different methodologies used; dNBR (differenced Normalized Burn Ratio) of short-term composites from Landsat images used alongside hotspot data was the most consistently successful in representing BA. We stress caution using GABAM in these regions, especially for the years 2001–2013, as Landsat-7 ETM+ scan lines are mistaken as burnt patches, increasing errors of commission. On the other hand, we highlight using regional products where possible, such as ABoVE-FED or ABBA in North America, and the Talucci et al. fire perimeter product in Russia, due to their detection of smaller fires which are often missed by global products.


Introduction
Located at high latitudes across North America, Fennoscandia, and Eurasia, the boreal forest, also known as the taiga, is the world's largest land biome, containing about 30% of the total forested area [1].The boreal forest is one of the world's largest biogeoclimatic areas, with an assortment of climates, topographies, fauna, and flora, making it one of the most important ecosystems on the planet [2,3].In North America, Fennoscandia, and Eurasia, the boreal forest is generally bounded by regions of Arctic tundra to the north and steppe and other temperate ecosystems to the south [4].
The treeless tundra landscape encompasses a range of dominant vegetation types including polar desert, tussock and shrub tundra, and wetlands [5].In many areas there may be permafrost present, either continuous, discontinuous, or sporadic, depending on climate, landscape position, vegetation, and disturbance history, among other factors [6].The steppe comprises mostly shrublands and grasslands [7].Both tundra and steppe vegetation result from prevailing limiting climatic conditions: low temperature in the tundra and water-energy balance relations in the steppe [8,9].
The Arctic has warmed 2-4 times faster than the rest of the planet since the 1970s [10,11], a phenomenon termed "arctic amplification".Consequently, boreal and tundra ecosystems are experiencing more frequent and higher intensity extreme weather events and disturbances, particularly wildfires [12][13][14].Burned area from wildfires and the associated greenhouse gas emissions have also increased during this time period [15,16].The huge extent of the Arctic-boreal zone means that any changes to the ecosystem, including those caused by increased emissions from methane, carbon dioxide, and aerosols, plus changes to land surface albedo, will have a significant impact upon the magnitude of future global climate change [17][18][19][20][21].
Fire has long been considered an integral part of the ecological process of boreal forests and remains the main disturbance today [22][23][24][25].Fire allows for large-scale germination and regeneration of component species, which may lead to significant, potentially longterm changes in forest biodiversity over time [26,27].Further, post-fire recovery of the taiga has been shown to be affected by both flora type and fire severity [28][29][30].Recently, intensifying fire regimes have been shown to have serious negative impacts through forest regeneration failure [25], and, at both regional and global scales, through the loss of native fauna, permafrost thaw, harm to human health, and damage to infrastructure [31][32][33].Although the majority of wildfires globally are caused by human activities [34,35], at more northerly latitudes with a low human population density, the majority of burned area in tundra and in boreal forests is caused naturally by lightning strikes [36][37][38].
Despite being mostly devoid of forested areas, wildfires still occur in regions of tundra and steppe, including overwintering fires, which also occur in taiga [39,40].Overwintering fires are those that ignite then smolder under the layer of peat and snow through the winter months, before re-emerging early the following spring [41].
Different land cover types give rise to different fire regimes when wildfires occur.For example, the large number of spruce and jack pine trees in North America promotes high intensity crown fires, whereas fires in Siberia are more often low-moderate intensity surface fires, due to greater coverage of larch trees, whose needles have a high moisture content, and Scots pines, which shed their lower branches and have thick bark, as well as more extensive shrublands and grasslands [42,43].
Wildfires in North America generally lead to high immediate tree mortality rates, with a considerable amount of debris falling to the ground, increasing the amount of fuel available to burn in future years [44].Surface fires in Siberia often cause extensive tree mortality that occurs over a series of years after the fire, with the changing climate, compounded by the timing of masting events in relation to time since fire, leading to delays in the regeneration of tree stands post-fire [29,45].Grass and shrubs often appear a few months after an area has been burned [46,47].Whether this is an initial stage of succession in a gradual return to forest dominance once again, or whether grasses and shrubs continue to prevail as a result of recruitment failure of tree species, is a key consideration as the climate warms and becomes drier [29].
Given the varied nature and significant impacts of wildfires across different ecosystems, satellite images are used to effectively detect and monitor global fire activity [48].Satellite images help determine the changing nature of our planet by identifying natural hazards such as landslides and floods, as well as forest fires [49][50][51].Images from three main groups of sensors, those onboard Landsat satellites, the Advanced Very-High Resolution Radiometer (AVHRR), and the Moderate Resolution Imaging Spectroradiometer (MODIS), have been used, alongside others such as SPOT and Sentinel-2, in the development of burned area (BA) products.
Products depicting BA from fires are generally developed from a combination of satellite imagery and active fire, or "hotspot", data, and derived land cover or vegetation maps.The hotspots are thermal anomalies, which can be caused by wildfires, smoke, agricultural fires, gas flares, or volcanoes [52].Each BA dataset is developed using differing algorithms at different spatial resolutions, which accounts for variations in the amount and location of BA, and scale and extent of the final product.The algorithms use various imagery input sources and consider different spectral indices.All, however, have to cope with the presence of clouds, which can complicate the process by obscuring the land for days, if not weeks, at a time and can be easily confused with smoke from fires [53,54].Full information on production of each of the BA datasets is given in the Methodology, Section 2.2.
Given the large number of sensors and resulting BA products, with different spatial and temporal resolutions, assessing BA products is critical for understanding changing fire regimes in northern ecosystems.While comparisons of BA products have been conducted on regional scales across different biomes, none have done so for the Arctic-boreal region.A list of previous relevant literature is given in Table 1, with a more complete review provided in the Discussion below.
Here we use 12 different BA datasets and MODIS active fire data to determine the total yearly BA and seasonal patterns of fires in Arctic-boreal North America, central/northern Siberia, and southern Siberia for the years 2001-2020.Based on this analysis, we also consider the advantages and limitations of each dataset.

Study Areas
Our study is focused upon three broad regions, one in North America and two in Russia, and is partly constrained by the data products we compare.In North America, the analysis was conducted over the entire boreal and tundra biomes in Alaska and Canada as defined by Olson et al. (2001) [62], whereas in Russia eight terrestrial ecoregions within the boreal and tundra biomes in the Siberian and Far Eastern Federal Districts were used, sensu Talucci et al. (2022) [63].Hereafter these regions are referred to as North America and central/northern Siberia, respectively.The eight ecoregions comprising central/northern Siberia are the Bering tundra, Cherskii-Kolyma Mountain tundra, Chukchi Peninsula tundra, East Siberian taiga, Northeast Siberian coastal tundra, Northeast Siberian taiga, Taimyr-Central Siberian tundra, and Trans-Baikal Bald Mountain tundra (see Talucci et al. (2022) [63] Figure 1).This area spans ca.7.75 million km 2 from ca. 50-78 • N and 79 • E-169 • W, whereas the North America area spans ca.8.86 million km 2 from ca. 46-83.5 • N and 172.5 • E-52.5 • W.
A second region was defined in Russia to investigate the performance of the global BA datasets in southern Siberia, which has a somewhat distinct fire regime when compared to central/northern Siberia [64,65].This latter region's northern boundary is the Southern Southern Siberia is the study region with the highest proportion of grasslands, known as steppe (19% of the total area), and mixed forests (29%) [68].Here there is no single dominant tree species, although there are several distinct regions of mixed forests.Birch, Scots pine, Siberian pine, and larch can be found throughout the whole region, with significant groups of Eurasian aspen (Populus tremula) in the west, and fir (Abies spp.) in the center [43].Vast areas of steppe are located between the forests throughout the entire southern Siberia region.Black spruce (Picea mariana) and jack pine (Pinus banksiana) are the two most dominant tree species that burn across North America [65,67].The former is located across both Alaska and Canada, whereas the latter exists only in boreal Canada.Water bodies and permanent wetlands are prominent across the region, especially in areas of Arctic tundra, comprising approximately 24% of the total area [68].Other major tree species across North America include various poplars (Populus spp.), balsam fir (Abies balsamea-Canada only), quaking or trembling aspen (Populus tremuloides), and white spruce (Picea glauca) [69].
In contrast to North America, larch (Larix spp.) is the dominant tree species in central/northern Siberia, with three species of pine (Siberian dwarf (Pinus pumila), Siberian (Pinus sibirica), and Scots (Pinus sylvestris)), plus birch (Betula spp.) constituting most of the remainder [43].The majority of fires (ca.82%) occur in larch and Scots pine forests [65].Larch tree cover spreads from the Yenisey river (the region's western border at ~90 • E) through the central part of Siberia to the far northeast of the country to ca. 165 • E. Siberian dwarf is the most common species of pine and is most prevalent in the eastern half of the region.Birch trees are found interspersed between areas of larch at mid and low latitudes.Scots pine is located at lower latitudes across the whole region, whereas there are two clusters of Siberian pines, in the south and southwest sections of the region, respectively.As with North America, there is an absence of trees in the tundra at the northern part of central/northern Siberia, although a large percentage of the region (ca.40%) is classified as open shrublands according to a global land cover product [68].
Southern Siberia is the study region with the highest proportion of grasslands, known as steppe (19% of the total area), and mixed forests (29%) [68].Here there is no single dominant tree species, although there are several distinct regions of mixed forests.Birch, Scots pine, Siberian pine, and larch can be found throughout the whole region, with significant groups of Eurasian aspen (Populus tremula) in the west, and fir (Abies spp.) in the center [43].Vast areas of steppe are located between the forests throughout the entire southern Siberia region.

Datasets
Twelve BA products of varying spatial and temporal resolutions, types, and extents were included in the analysis, along with MODIS active fires.A list of the products used can be found in Table 2, followed by summaries of the methodologies of each dataset.One method which appears frequently is the differenced Normalized Burn Ratio (dNBR), which is the change between the pre-and post-fire Normalized Burn Ratio (NBR).NBR is calculated using the Near Infrared (NIR) and Short-Wave Infrared (SWIR) bands from satellite images using Equation (1).[84]) and FireCCI50 (2001( -2016( : Chuvieco et al. (2018) ) [85]).The product uses data from four different input sources to generate BA information.The Red and NIR bands were taken from the MODIS Terra MOD09GQ (Collection 6) product [86] (250 m spatial resolution), with additional pixel state data from the Terra MOD09GA product (1 km).The 1 km MODIS Global Monthly Fire Location Product, MCD14ML (collection 6), was used to extract active fire (hotspot) information.Finally, the 300 m Land Cover CCI v2.0.7 product [87] was used to generate a burnable mask for each year, which was created by removing urban and bare areas, water bodies, and permanent snow and ice classes.Subsequently, monthly image composites were made from cloudless and cloud shadow-less images, based on the date of hotspot detection and the minimum NIR value.The minimum NIR was taken as the lowest value within a 10-day window preand post-fire, with up to 15 post-fire days included if less than four valid observations were obtained.Hotspot spatial clusters were then generated, with a radius of influence of 1875 m (1500 m buffer with 375 m to account for potential geolocation issues), and a temporal threshold of 4 days.These clusters formed the basis of the seed phase, where only pixels with a high probability of being burned were selected.A growing phase was then carried out on these pixels under certain threshold conditions.Finally, patch filtering was conducted to reduce commission errors caused by excessive growing and remove noise before the final product was created.FireCCI51 reports a commission error (false positives included) of 54.4% and an omission error (false negatives included) of 67.1% for the years 2003-2014.
The FireCCI51 product is available at two spatial resolutions projected over the WGS84 reference ellipsoid: pixel (250 m spatial resolution) and grid (0.25 • ), with monthly global BA data from 2001 to present.The grid product is formed by summing the number of burned 250 m pixels within each 0.25 • grid.For this analysis we chose to use the pixel product, as the day of burn and a confidence score are also given for each pixel in separate bands.We initially considered all pixels with a 50% or greater confidence of being burned; however, after preliminary analysis, this was subsequently increased to greater than or equal to 80% (Figure S1) to reduce the likelihood of commission errors.The value of 80% was chosen to align with the upper confidence bound of the MODIS active fires.

MCD64A1
The Collection 6 (C6) MCD64A1 product [71] is the third 500 m BA product released by NASA, following the Collection 5.1 (C5) MCD45A1 and MCD64A1 products.MCD64A1 takes advantage of the C6 suite of MODIS products, including bands 1, 5, and 7 from the Terra/Aqua 500 m atmospherically corrected Level 2G daily surface reflectance products (MOD09GHK/MYD09GHK, respectively [88]), 1 km Terra/Aqua level-3 daily active fire products replicated to 500 m resolution (MOD14A1/MYD14A1 [89]), and the 500 m Level-3 annual land cover product MCD12Q1 [90].First, the surface reflectance time series values were extracted such that the observations were cloud-free, fire-free, over land, and had physically valid reflectance.If multiple values for each day were available, then the observation sensed at the smallest view zenith angle was chosen.Grid cells having too few valid observations to fill the pre-and post-burn windows were labeled as unclassified and subsequently ignored.Subsequently, active fire composite imagery was produced based on the maximum separability of pre-and post-fire burn-sensitive vegetation index (VI) values over sliding 8-day pre-and post-fire windows.Several burned and unburned grid cells were then selected to be used as part of the model training process to estimate probabilistic thresholds suitable for classifying individual 500 m grid cells as either burned or unburned.The probabilities for each cell of being burned prior to or after a fire event were estimated, given the observed change in VI, before being classified subject to the previously determined thresholds and labeled.MCD64A1 reports a commission error of 24% and an omission error of 37%.
The C6 MCD64A1 monthly BA product is available globally from 2001 to present at a 500 m resolution on the MODIS sinusoidal projection, with the approximate day of burn given for each pixel.A second "Uncertainty" band was added which indicates the uncertainty of the burn date (±number of days); however, this was not used in this analysis.

GABAM
The Global Annual Burned Area Maps product (GABAM) [72] used 6 bands (Blue, Green, Red, NIR and two SWIR bands) across all available Landsat-8 Surface Reflectance Tier 1 and Tier 2 collections as its core data source.Pixels masked by clouds, cloud shadows, water, snow, ice, or filled/dropped pixels were excluded, leaving only clear land pixels.Land cover information was extracted from the MODIS MCD12C1 product [91], from which 7 land cover categories were ultimately considered.The BA density of the GFED4 product [73] was used to create five equal-frequency fire intervals with quantile classification.Training and validation sites were subsequently chosen by stratified random sampling with a 3:2 ratio.Alongside the surface reflectance of the 6 bands, 8 spectral indices (2 NBR indices, Burned Area Index, Mid-Infrared Burn Index, Normalized Difference Vegetation Index (NDVI), Global Environmental Monitoring Index, Soil-Adjusted Vegetation Index and Normalized Difference Moisture Index) were also derived to perform the global BA mapping.A random forest (RF) algorithm was used to train a decision forest classifier based on the aforementioned input features.The outputs of the RF algorithm were used alongside NDVI and dNBR thresholds developed from imagery of the current and previous year to select seeds to be used as the starting point for growing each fire.Finally, an iterative procedure around each seed pixel was used to grow each fire based on the burned probabilities of the eight neighboring pixels being greater than or equal to 0.5 to give the final BA product.GABAM reports a commission error of 13.2% and an omission error of 30.1%.
Note that the methodology for the GABAM publication is given specifically for the year 2015.The collection of Landsat-8 images was used as its principal data source; however, Landsat-8 was only launched in February 2013, whereas the GABAM product covers the entire 1985-2021 period on an annual temporal resolution.Landsat-5 and Landsat-7 images were used prior to 2013 to ensure a consistent 30 m BA product; however, there is no mention of this in the paper, nor if the authors corrected for the Landsat-7 Scan Line Corrector failure in May 2003.

GFED4s
The Global Fire Emissions Database version 4.1 with small fires (GFED4s) is an extension of the original GFED4 product [73], with the additional "s" signifying "small fires".Both GFED4 and GFED4s also give information about carbon emissions from fires alongside BA [74].The methodology to identify small (sub-500 m) fires was first given in Randerson et al. (2012) [92].GFED4s begins in 1997 and ends in 2016; therefore, two different methodologies were used depending on whether MODIS data were available (August 2000 onwards) or not.During the MODIS era, the 500 m MCD64A1 Collection 5.1 product, aggregated to a 0.25 • grid, was used as a baseline.This product was then combined with 1 km MODIS Terra/Aqua thermal anomaly data and 500 m surface reflectance observations to estimate additional BA associated with small fires in each grid cell.The MODIS fire location product MCD14ML was additionally used to improve geolocation accuracies, and only detections with small or moderate scan angles were retained to reduce uncertainties in the geolocation.dNBR was then estimated for those active fire detections for which all four overlapping pixels were also classed as burned.Note that data from the Aqua satellite were only available from January 2003, so for the period August 2000-December 2002 a climatological ratio was used to estimate the small fires alongside data from Terra.In the pre-MODIS era (January 1997-July 2000) active fire data were extracted from either the Visible and Infrared Scanner (VIRS) or the Along Track Scanning Radiometers (ATSRs) to estimate BA.VIRS only came online in January 1998, so for the first year only ATSR data were used.A relationship was developed between aggregated active fires in the pre-MODIS era and the BA during the MODIS era using linear regression models, before the aggregated BA was distributed to individual 0.25 • grid cells.The final GFED4s BA product is global with a monthly temporal resolution and 0.25 • spatial resolution.

GFED5
GFED5 [75] takes advantage of the availability of high spatial resolution satellite data from Landsat (30 m) and Sentinel-2 (20 m), as well as the updated Collection 6 MODIS suite, including the MCD64A1 BA product, 500 m MCD12Q1 land cover product [93], and 250 m Terra and Aqua combined Vegetation Continuous Fields product MOD44B [94].As with GFED4s, the Terra/Aqua MODIS active fire products were used from 2001, with ATSR/VIRS data again utilized in the pre-MODIS (1997-2000) era.MCD64A1 C6 was used as the core BA input during the MODIS era, with nine Landsat or Sentinel-2 datasets used for calibration and an additional eight for validation.To account for commission errors of unburned "islands" within MCD64A1 pixels, BA was multiplied by a correction scalar.Similarly, omission errors caused by the absence of small fires or incompletely mapped burned areas within MCD64A1 were accounted for by multiplying a second corrective scalar with the surface area of MODIS active fire detections outside the perimeter of the 500 m MCD64A1 pixels.These scalars were derived from the fine-resolution images and the MODIS BA and active fire data.Note that the BAs for croplands, tropical peatlands, and deforestation regions were calculated separately before being summed together with the normal land cover BA.The final GFED5 global BA product for 2001-2020 was produced at a 0.25 • spatial resolution and monthly temporal resolution.
The methodology for the pre-MODIS era is similar to that of GFED4s, although the data are produced at 1 • , as opposed to 0.25 • , due to poor spatial variability of BA derived from the ATSR and VIRS active fires.

FireCCILT11
This product was derived by the ESA from the Land Long Term Data Record (LTDR), which provides a continuous dataset of AVHRR images that have been corrected both geometrically and radiometrically [76].The imagery was collected by the AVHRR 2/3 sensors aboard NOAA-7, 9, 11, 14, 16, 18, and 19 satellites before being stored at a reduced 4 km resolution.The LTDR resamples these data to 0.05 • (approximately 5.5 km) and the Surface Reflectance Product of this dataset was used.Note that although the final product is provided from 1982 to 2018, the year 1994 was removed due to an absence of several dates in the LTDR dataset; 2019 and 2020 were also excluded as a result of poor quality LTDR data.Land cover data from the Copernicus Climate Change Service [95] and the ESA Climate Change Initiative [87] were used to create a burnable mask for each year from 1992, the same process as with FireCCI51.For the pre-1992 period, the land cover for the year 1992 was used.Monthly imagery compositing was employed to mitigate problems from daily image acquisitions and cloud-free images were removed.An LTDR burned index (LBI) was developed based on the input channels and derived bands.The bands considered were Red, NIR, medium infrared reflectance, brightness temperature for channels 3, 4, and 5, view and solar zenith angles, relative azimuth, and a quality assessment layer.Two different random forest models, one global model and a second specifically for the boreal region, were developed to detect BA and trained using FireCCI51 for the years 2001-2018.The boreal model was developed as a result of poor estimations in boreal regions of the global model.The global model considered composites using the LBI 5 months pre-fire and 6 months post-fire, whereas the window for the boreal model was reduced to 3 months before and 4 months after fire occurrence.Detection thresholds were used to create binary burned/unburned pixels before an estimation of the total BA in each burned pixel was performed for both models.The final BA product was achieved by combining the global and boreal models into one unified product.
Both FireCCILT11 and FireCCI51 are produced by the ESA and as such share several formatting similarities.As with FireCCI51, FireCCILT11 is available as both a pixel (0.05 • ) and grid (0.25 • ) product, and the pixel dataset offers an uncertainty band.We again initially considered all pixels with a confidence of greater than or equal to 50%; however, this was subsequently decreased to 30% or greater due to poor performance of the BA product (Figure S2).The 30% confidence value was chosen to align with the lower confidence bound of the MODIS active fires.For this analysis we used the pixel product as the approximate day of burn is provided for each pixel.

MTBS
The first region-specific and fire polygon-based product considered is Monitoring Trends in Burn Severity (MTBS) [77], covering the United States of America at a 30 m spatial resolution and generated annually from 1984 to present.Imagery from the Landsat 4, 5, 7, and 8 satellites was used as the core input, and Landsat 9 images will be considered from 2023.The NBR of pre-and post-fire images was calculated using the NIR and SWIR bands and dNBR was computed.The dNBR images were then used to determine the severity of each fire and to delineate the fire perimeters.This last step relied on the interpretation of several analysts using their past experience of fire behaviour and effects.The final perimeter product also includes the start date for each polygon.

NBAC
The National Burned Area Composite (NBAC) [78] is compiled using three different BA products.The first is from the Canadian National Fire Database (CNFDB), where fire polygons are mapped annually by several fire management agencies.However, due to inter-annual variation caused by a lack of a unified mapping system, the NBAC considers these polygons only as potential sources which may or may not be used in the final product.
The second product is the Multi-Acquisition Fire Mapping System (MAFiMS), using atmospherically corrected Landsat 5, 7, and 8 pre-(1 year) and post-burn imagery.Pixels masked by cloud or covered by water were excluded, as were agricultural activities.The Red and NIR bands were used to calculate NDVI, which was compared before and after fire occurrence to define unburned regions.dNBR was then determined for each image to produce a threshold, which identified regions of high vegetation change, where some pixels will likely have been burned.This high-change layer was then converted to polygons, buffered at 1500 m, and filtered using MODIS hotspot data to create units which fully contain individual fire events.Active fire data were further used to determine the dates in which each event occurred to ensure they took place only in the year of interest.Finally, for quality control, each polygon was visually inspected and adjusted if necessary.
The third data source used to determine NBAC is the Hotspot and NDVI Differencing Synergy (HANDS) algorithm [96], which provides polygons for large (>250 ha) fires not reported either by the agencies or through using MAFiMS.The HANDS algorithm used 1 km SPOT VEGETATION imagery, which was replaced by Prova-V imagery resampled to 250 m in 2013.The algorithm used a normalized NIR-SWIR index based on pre-and post-burn imagery alongside active fire and land cover data.Burned pixels were then mapped based on given thresholds and were converted to polygons to be incorporated into NBAC.As coarse satellite images tend to overestimate the total BA [97], a Theil-Sen regression [98] was used alongside area-based calibration models and spatial and temporal filters to generate statistically unbiased estimates of the BA from HANDS polygons.
The methodology described above was taken from Hall et al. (2020) [78], referring to fires for the period 2004-2016.This was later extended in Skakun et al. (2022) [79] to cover the years 1986-2020 and is summarized here.CNFDB polygons for the early period were obtained from the Canadian Wildland Fire Information System [99] using various methods, including manual sketch mapping and aerial GPS, which are often far less accurate than Landsat-derived polygons [100].MAFiMS was extended in recent years using the same methodology but also incorporated Sentinel-2 satellite imagery.MAFiMS was extended back using a combination of AVHRR and MODIS active fire detections and Landsat collection 1 archive imagery.Finally, polygons from the HANDS algorithm could be developed back to the start of SPOT VEGETATION satellite data in 1995.A fourth product, the National Terrestrial Ecosystem Monitoring System Composite to Change [101], was also integrated to better represent some fire events that were missed by MAFiMS.
The complete NBAC polygon product is continually updated annually at a 30 m spatial resolution across Canada.The ignition date of each fire polygon occurrence is also given.Note for this analysis the MTBS product for boreal/tundra Alaska is used alongside NBAC for boreal/tundra Canada to create a combined MTBS/NBAC polygon product.

ABoVE-FED
The Fire Emissions Database (FED) from NASA's Arctic Boreal Vulnerability Experiment (ABoVE) program contains BA information for the entirety of Alaska and Canada from 2001 to 2019 at a 500 m spatial resolution and annual timescale [80].Atmospherically corrected cloud-free Landsat 5, 7, and 8 30 m images were extracted and 1-year pre-and post-fire NBR values were calculated at field-based burned and unburned control sites to create a dNBR threshold.A 30 m resolution unburnable land cover mask for each year was created using the Joint Research Center's yearly water history product [102] and the 2010 land cover product of the North American Land Change Monitoring System [103].A 30 m binary burned/unburned Landsat product was generated using the dNBR threshold and vegetation mask.The product was then upscaled to the MODIS 500 m resolution and projection, assigning the value of the percentage of burnable area within each 500 m pixel as the associated burn fraction when >50% of the 30 m vegetated pixels were burned within it.This process took place when at least 85% of the Landsat data were available, accounting for 81% of the total burned pixels.For the remaining 19%, the pixel was classified using a similar process derived from MODIS Terra (MOD09GA) and Aqua (MYD09GA) imagery.To account for fire detection issues arising from using imagery 1-year pre-and post-fire, a seasonal MODIS product was also developed with a 30-day pre-and post-fire window, masking out pixels with 10 or fewer valid observations in each window.The day of burn for each pixel was taken from MODIS active fire data where possible, with an interpolation approach or the MODIS seasonal BA product used where no active fire data were associated with a given polygon.ABoVE-FED reports a balanced commission/omission error of 6.6% when using Landsat data and 14.2% using MODIS data.

SRBA
The Surface Reflectance Burned Area (SRBA) product [81] is derived from two MODIS products: MOD09 daily surface reflectance data and MOD14 active fires.Pre-processing was carried out at the MODIS 500 m spatial resolution before being resampled to 250 m.Pixels with a high view or solar zenith angle were first masked out, and then pixels contaminated by clouds were detected using threshold values from the Normalized Difference Snow Index, using the Blue and SWIR bands from the MOD09 product, and also masked out.NIR was analyzed along the cloud-shadow line to correctly identify areas covered by cloud shadow, which were then removed.A moving monthly time window was used to filter out residual noise.dNBR was calculated based on daily composites from uncontaminated pixels, comparing the NBR in pre-and post-fire seasons.Note that although Bartalev et al. (2012) [81] call this variable the Normalized Short-Wave Vegetation Index, the formula used is identical to NBR.Potentially burned pixels were identified using a dNBR threshold, and spatial filtering was used to remove pixels affected by non-fire disturbance factors.Finally, pixels were grouped together into polygons for each day by comparing fire events with MODIS active fire data.
SRBA is produced at a daily temporal resolution covering the entirety of the Russian Federation from 2006 to present.We were provided with access to SRBA data for longitudes from 60-180 • E for the years 2006-2020.The data were given to us with an annual temporal resolution, so the day of burn for each pixel was not provided.

Talucci et al. Fire Perimeter Product
The Talucci et al. (2022) [63] fire perimeter product covers eight boreal and tundra ecosystems in Russia, as outlined in Section 2.1.Landsat 5, 7 and 8 Surface Reflectance Tier 1 images during the growing season months (May to October) were chosen to ensure a high spatial resolution (30 m).Pixels covered by snow and clouds were masked, before annual composite images were generated using the best pixel approach [104] to allow for missing data from snow and cloud cover, and the Landsat-7 scan line error.dNBR was calculated using one-year pre-and post-fire composites before each image was binary classified as either burned or unburned.MODIS active fire data were used to create polygons, before the fire perimeters were delineated around the corresponding Landsat pixels and exported as shapefiles.The polygons were buffered at a 1 km resolution to combine adjacent polygons before removing any holes and applying the reverse buffer.A smoothing and simplification function was then applied to reduce the pixelated appearance and the number of vertices of each polygon.MODIS active fires were again used to estimate the polygon start and end date and to eliminate perimeters with a 0% confidence rating.Small fires (<200 ha) were removed to match the criteria in the Canadian large fire database [24] but low-confidence fires were retained as Siberian fires often burn at a lower intensity.
The final product was produced at a 30 m resolution from 2001 to 2020.When considering the seasonality of fires in our analysis, we used the MinDay value (the given start date for each polygon).The difference between MinDay and AvgDay (the mean date in between the start and end dates) is shown in Figure S3.

ABBA
The Arctic Boreal Burned Area (ABBA) [82] dataset is provided at a 500 m resolution covering the boreal and tundra biomes on an annual timescale for the years 2002-2022 and is designed to capture fires occurring later in the season.Several MODIS products were used as inputs, including the C6 MCD14ML active fires [83], C6 MOD09A1 Surface Reflectance [105], C5 MOD44W land-water mask [106], and the C5 MOD44B Vegetation Continuous Fields Collection [107].Eight-day surface reflectance composites from MOD09A1 were collected in the year of the fire alongside that of the previous year.Poorquality pixels were masked out and dNBR grids were produced best on an 8-day preand post-fire window.Thresholds were then developed based upon vegetation indices to separate the dNBR grids into binary "potentially burned" and "unburned" categories, before being compared with the active fire data to eliminate any changes due to non-fire reasons.The resultant 8-day burn masks were used to create the final BA product.

MODIS Active Fires
Although MODIS active fires are used in this analysis, they should not be regarded as a substitute for a BA dataset, as the relationships between active fires and burned areas are not spatially or temporally constant [108,109].However, the active fire data can help bound BA estimates and are particularly useful for diagnosing seasonality.Here we use MOD in this summary to refer to both Terra (MOD) and Aqua (MYD) as the detection algorithm remains the same for both satellites.Note that we chose to use active fire data from MODIS rather than VIIRS-not to be confused with VIRS-as MODIS data are provided for the whole time period, whereas VIIRS was only launched in October 2011.
Several bands of different wavelengths from MOD021KM Swath Radiance and MOD03 Swath Geolocation data from the satellites were extracted, aggregated to a 1 km spatial resolution, checked and used as inputs into the algorithm.Coastal pixels were detected and removed from consideration and images obscured by clouds were also removed.Daytime fire pixels were tentatively identified using thresholds derived from the input data; however, the final identification for nighttime fire pixels was conducted at this stage.
Further rejection tests were carried out for water and daytime land pixels based on more thresholds, and include sun glint rejection, desert boundary rejection, land-pixel coastal rejection, forest clearing rejection, and water-pixel coastal rejection.Finally, each pixel which passed these tests was assigned a confidence rating.In our analysis we considered pixels with both a nominal (30-80%) and high (>80%) confidence to include the low intensity fires in southern Siberia.Note there was a discontinuity in July 2002 when Aqua came online, substantially increasing the total BA.There were also major science data outages for both Terra (15 June-3 July 2001), due to a MODIS Characterization Support Team outage, and Aqua (16 August-2 September 2020), caused by a Formatter-Multiplexer Unit/Solid State Recorder error.

Analysis Techniques 2.3.1. Yearly BA
The total BA for each year during the MODIS era (2001-2020) was computed using the sum reducer (ee.Reducer.sum[110]) in Google Earth Engine (GEE) [111].The datasets were either uploaded as assets into GEE or imported directly from the Earth Engine Data Catalog.Each pixel-based product was uploaded as a collection of GeoTIFF images and therefore processed as an ee.ImageCollection, whereas each polygon-based product became an ee.FeatureCollection as the datasets were uploaded as shapefiles.
For those datasets where the whole pixel is designated as burned (FireCCI51, MCD64A1, GABAM, SRBA, ABBA, MODIS active fires), the burned pixels were given the value 1, with unburned pixels assigned 0. The sum reducer then summed all the pixels of value 1 in a given region, which was then multiplied by the area of the pixel to give the total BA in square meters (m 2 ), before being converted to megahectares (Mha).
A similar technique was applied to those datasets which provide the burned fraction of each pixel (FireCCILT11, GFED4s, ABoVE-FED).Here the original value of each pixel was maintained before applying the sum reducer, multiplying by the area of the pixel and converting to Mha in each region.
The total BA in GFED5 was computed using a variation of this process: in the dataset each pixel was given the value of the BA in square kilometers (km 2 ), so the sum reducer was applied before being converted to Mha.
As the polygon-containing datasets (MTBS, NBAC, Talucci et al.) were uploaded as shapefiles to GEE, they were necessarily processed differently to the pixel-based products.The area of each polygon in m 2 was computed using the geometry.areafunction in GEE which was mapped over the whole region, before the total area was summed over each polygon using the sum reducer and subsequently converted to Mha.
This study used two methods to counteract the possibility of repeated fire events taking place in a single year.For pixel datasets where the approximate day of burn was provided (FireCCI51, MCD64A1, FireCCILT11, and ABoVE-FED), we assumed that each pixel could only be burned once in a given year, based upon the fire return intervals in the Arctic-boreal region [42,112].Alternatively, for GFED4s, GFED5, and MODIS active fires we conducted our analysis monthly to ensure that repeat events could occur on the same grid cells or pixels within a given year, and all would be accounted for.Polygon products and pixel datasets with an annual temporal resolution were not affected by repeat fire events.

Monthly Fire Percentage
The monthly proportion of each year's total fires (e.g., in 2004 around 45% of the fires in North America occurred in July) was also computed in GEE using the unweighted frequency histogram reducer (ee.Reducer.frequencyHistogram.unweighted [113]).This could only be computed where either the temporal resolution of the dataset was monthly or the Julian Day value was provided (Table 2).
For the pixel datasets with a daily temporal resolution (FireCCI51, MCD64A1, ABoVE-FED, FireCCILT11, MODIS active fires), along with Talucci et al., the first day in the year of fire occurrence detection, often referred to as the Julian day, was given as a number between 1 and 365 (366 for leap years).The unweighted frequency histogram reducer was then used to find the total number of pixels with the same Julian day, and subsequently the total daily BA, before grouping by month.The proportion of each month's fires was then calculated as a percentage relative to the whole year.The MTBS and NBAC data were processed similarly, although these had already been grouped by month, bypassing the first step.
For GFED4s, GFED5, and MODIS active fire data, the total BA was computed monthly, then divided by the yearly total to find the relative percentage.
When referring to seasonal patterns, we split the seasons as follows: winter (December, January, February), spring (March, April, May), summer (June, July, August) and autumn (September, October, November).The method for calculating the percentage difference between datasets is given in Appendix A.

North America
The BA for each dataset per year from 2001 to 2020 is given in Table S1 and depicted in Figure S4.However, the analysis for North America focusses upon the years 2002-2016 in order to include all the datasets.The total BA for each dataset during the analysis period is given in Figure 2a, with the related box and whisker plots in Figure 2b.The largest fire year in North America during this period was 2004, with an average BA of 4.88 ± 0.89 Mha (standard deviation (sd)) across the products, with a mean BA across the datasets of 2.55 ± 1.24 Mha/year.On average, more area is burned in the second half of the time period, 2.81 ± 1.34 Mha/year for 2010-2016, than the first, 2.32 ± 1.19 Mha/year for 2002-2009, demonstrating an increasing tendency.2004 was also the largest fire year over the entire 20-year period.
MODIS active fires are perhaps expected to produce the greatest amount of BA as the MODIS sensors detect all thermal anomalies, which can be related to non-wildfire activities such as cropland fires, partial burns, or even increased haziness across satellite images [114,115].However, in North America this is not the case.MODIS active fires detect significantly less BA in Alaska and produce the greatest amount of BA in only 8 of the 20 years in Canada.Here most fires are high intensity, rapidly spreading crown fires, and as such MODIS sensors can miss burned pixels.MODIS active fires have also been known to struggle in this region due to fires occurring on cloudy days going undetected, underestimating the total number of fires by at least 18% [116].
Note that the datasets cannot be split into any statistically significant groups based upon average annual BA from 2002-2016 (Figure 2b and Table S2), as described in Appendix B. GABAM is the product which shows the least amount of BA between 2002 and 2016, with a total BA of 27.9 Mha.Datasets with a high spatial resolution, such as GABAM, will often produce less BA than those with lower spatial resolutions, as they will more precisely identify the edges of fires and unburned "islands" within them.However, the combined MTBS/NBAC polygon product, also at 30 m resolution, contains 46.8% more BA (Figure S5), meaning that differences in total BA are also related to the different product methodologies.The two pixel datasets developed specifically for North America, ABBA and ABoVE-FED, show the highest amount of BA, 46.6 and 48.3 Mha, respectively.In North America, six of the eight datasets where the day of burn was provided agree on the seasonal patterns of boreal fires for the period 2001-2016 (Figure 3), exhibiting a single peak in July with a three-month increase from April to July and a three-month decrease from July to October.All the datasets agree that the peak occurs in July, with ca.41% of all BA occurring during this month.With the exception of FireCCI51 and FireCCILT11, the datasets concur that June is the second biggest fire month, when an average of 33% of the yearly BA occurs, followed by August with ca.15%.Among these datasets there is more variation in June (10.5%)than August (4.4%)due to the difference in methodology between the pixel and polygon-based datasets.The MTBS/NBAC combined datasets show a tendency towards earlier occurrence in June followed by a steeper drop in August as the Julian Day value is taken as the first day of fire occurrence in each polygon, whereas the pixel datasets track the constant progression of the fire on a pixel-by-pixel basis, with later Julian Day values as each fire spreads.For FireCCI51, the proportions of June/August fires are decreased/increased, respectively, with an almost identical percentage, ca.22%, of BA occurring in both early and late summer.The FireCCILT11 dataset agrees that most BA in this region is produced in July (42%).However, the seasonal pattern otherwise differs significantly from the other datasets.FireCCILT11 shows that BA constantly increases by 2.5-3% per month from March to June, before jumping to the July peak, a small drop to August (37%), followed by a second, steeper drop down to 4% in September, suggesting that BA is produced toward the end of summer rather than the start.

Central/Northern Siberia
The analysis for central/northern Siberia spans 2006-2016 to include all datasets; however, the results for the entire 2001-2020 period are given in Table S3 and shown in Figure S6.The total BA for each dataset during the shorter analysis period is given in Figure 4a, with the related box and whisker plots in Figure 4b.Despite this study region having a smaller area than the North America region, an average of 4.04 ± 1.49 (sd) Mha/year burned in central/northern Siberia between 2006 and 2016, 58% higher than North America.The largest fire year during this period was 2014, with a mean BA across the datasets of 5.92 ± 2.17 Mha; however, over the full 20-year period it only ranks as the fifth largest after 2003 (14.For FireCCI51, the proportions of June/August fires are decreased/increased, respectively, with an almost identical percentage, ca.22%, of BA occurring in both early and late summer.The FireCCILT11 dataset agrees that most BA in this region is produced in July (42%).However, the seasonal pattern otherwise differs significantly from the other datasets.FireCCILT11 shows that BA constantly increases by 2.5-3% per month from March to June, before jumping to the July peak, a small drop to August (37%), followed by a second, steeper drop down to 4% in September, suggesting that BA is produced toward the end of summer rather than the start.

Central/Northern Siberia
The analysis for central/northern Siberia spans 2006-2016 to include all datasets; however, the results for the entire 2001-2020 period are given in Table S3 and shown in Figure S6.The total BA for each dataset during the shorter analysis period is given in Figure 4a, with the related box and whisker plots in Figure 4b.Despite this study region having a smaller area than the North America region, an average of 4.04 ± 1.49 (sd) Mha/year burned in central/northern Siberia between 2006 and 2016, 58% higher than North America.The largest fire year during this period was 2014, with a mean BA across the datasets of 5.92 ± 2.17 Mha; however, over the full 20-year period it only ranks as the fifth largest after 2003 (14.29 ± 5.31 Mha), 2020 (9.61 ± 4.00 Mha), 2019 (7.98 ± 3.19 Mha), and 2002 (7.44 ± 3.66 Mha).There is a large increase in the mean BA per year in the second half of the analysis period (2012-2016) compared to the first half (2006-2011), with fires increasing from 3.21 ± 1.29 Mha/year to 5.03 ± 1.10 Mha/year, and 4.39 ± 3.95 Mha/year to 5.64 ± 1.92 Mha/year for the 2001-2010 and 2011-2020 time periods, respectively.
In central/northern Siberia there is considerably more disagreement between the datasets compared to North America (Figure S7) and, unlike in North America, we cannot rely upon official statistics from government agencies, which have consistently been grossly underestimated [45,117,118].One primary reason for this is that boreal forest areas in the Siberian and Far Eastern Federal Districts, comprising 21.6% of the total forest area in Russia, are classed as "remote unused forests" and are excluded from official statistics [119][120][121].In central/northern Siberia there is considerably more disagreement between the datasets compared to North America (Figure S7) and, unlike in North America, we cannot rely upon official statistics from government agencies, which have consistently been grossly underestimated [45,117,118].One primary reason for this is that boreal forest areas in the Siberian and Far Eastern Federal Districts, comprising 21.6% of the total forest area in Russia, are classed as "remote unused forests" and are excluded from official statistics [119][120][121].We split the datasets into three distinct groups based on the computed BA statistical significance values (Figure 4b and Table S4), as described in Appendix B. Once again GABAM is the dataset with the lowest amount of average annual BA, 1.79 ± 0.77 Mha/year from 2006 to 2016, whereas MODIS active fires shows the most amount of average BA in this period: 6.80 ± 2.58 Mha/year.The group with the most amount of BA (Talucci et al., GFED5 and MODIS active fires) has over 3.5 times more average annual BA than GABAM.However, large differences can occur between individual datasets each year.For example, Additionally, we also found significant issues with GABAM during the years 2001-2013.Pixels are often incorrectly identified as burned in the form of additional "stripes" and/or "boxes" (Figure 5a), or scan lines appear as missing patches through burned areas (Figure 5b).The stripes are a result of Landsat-7 ETM+ scan lines being mistaken as burnt patches, increasing commission errors [122].The Landsat-7 ETM+ Scan Line Corrector (SLC) failed permanently on 31 May 2003, causing all subsequent images to miss ca.22% of the normal scene area [123].However, we found that boxes can occur if the GABAM processing algorithm fails at the point of overlap between two or more Landsat-7 ETM+ images.While we found that although these issues can occur all over the world, they are most prevalent in Siberia.
significance values (Figure 4b and Table S4), as described in Appendix B. Once again GABAM is the dataset with the lowest amount of average annual BA, 1.79 ± 0.77 Mha/year from 2006 to 2016, whereas MODIS active fires shows the most amount of average BA in this period: 6.80 ± 2.58 Mha/year.The group with the most amount of BA (Talucci et al., GFED5 and MODIS active fires) has over 3.5 times more average annual BA than GABAM.However, large differences can occur between individual datasets each year.For example, in 2007 Talucci et al. has 12.5 times more BA than SRBA (5.56 Mha vs. 0.44 Mha), and in 2012 ABBA has 11.8 times more BA than GFED4s (9.01 Mha vs. 0.77 Mha).
Additionally, we also found significant issues with GABAM during the years 2001-2013.Pixels are often incorrectly identified as burned in the form of additional "stripes" and/or "boxes" (Figure 5a), or scan lines appear as missing patches through burned areas (Figure 5b).The stripes are a result of Landsat-7 ETM+ scan lines being mistaken as burnt patches, increasing commission errors [122].The Landsat-7 ETM+ Scan Line Corrector (SLC) failed permanently on 31 May 2003, causing all subsequent images to miss ca.22% of the normal scene area [123].However, we found that boxes can occur if the GABAM processing algorithm fails at the point of overlap between two or more Landsat-7 ETM+ images.While we found that although these issues can occur all over the world, they are most prevalent in Siberia.6, the datasets show further disagreement in central/northern Siberia compared to North America.GFED4s, GFED5, and MODIS active fires show a constant increase to an initial peak (ca.16%) in May, followed by a drop to 12.5% in June, a second, larger peak ranging from 28 to 38% in July, and a steady decrease through August to ca. 5.5% in September, before dropping below 1.5% in October.When considering the seasonal patterns over the 2001-2016 period, as shown in Figure 6, the datasets show further disagreement in central/northern Siberia compared to North America.GFED4s, GFED5, and MODIS active fires show a constant increase to an initial peak (ca.16%) in May, followed by a drop to 12.5% in June, a second, larger peak ranging from 28 to 38% in July, and a steady decrease through August to ca. 5.5% in September, before dropping below 1.5% in October.
FireCCI51 and FireCCILT11 display consistent seasonal tendencies with each other but, as in North America, show an inflated proportion of late summer BA compared to early summer.These products again have an initial small peak in May (ca.17.5%); however, the main summer peak comes in August (ca.30.7%), rather than July (ca.25.4%), before dropping sharply to 5.8% in September and decreasing to the end of the year.
MCD64A1 is anomalous for the pixel-based products in that it does not recognize an initial spring peak, with BA spread evenly between March and June; approximately 8.5% of the yearly BA happens in each of the four months.This is then followed by a jump to the July peak at 37.1%, before dropping through 20.8% in August to 7.3% in September and a negligible amount for the rest of the year.FireCCI51 and FireCCILT11 display consistent seasonal tendencies with each other but, as in North America, show an inflated proportion of late summer BA compared to early summer.These products again have an initial small peak in May (ca.17.5%); however, the main summer peak comes in August (ca.30.7%), rather than July (ca.25.4%), before dropping sharply to 5.8% in September and decreasing to the end of the year.
MCD64A1 is anomalous for the pixel-based products in that it does not recognize an initial spring peak, with BA spread evenly between March and June; approximately 8.5% of the yearly BA happens in each of the four months.This is then followed by a jump to the July peak at 37.1%, before dropping through 20.8% in August to 7.3% in September and a negligible amount for the rest of the year.
The seasonality of the Talucci et al. polygon product can be compared to that of the MTBS/NBAC polygon product in North America, as it shows a tendency towards earlier fire occurrence as it does not take into account the spread of the fire between pixels, rather giving the date of the initial fire for the whole polygon.This can be somewhat improved by using the AvgDay value, which accounts for an increased proportion of late summer BA.However, both products show a single peak at approximately 35.4% in July following a constantly increasing number of fires over the previous four months and preceding a rapid two-month drop.

Southern Siberia
For southern Siberia the analysis is limited to the six global BA datasets (FireCCI51, MCD64A1, GABAM, GFED4s, GFED5, FireCCILT11) and MODIS active fires for the years 2001-2016.The results for the whole 2001-2020 period are shown in Table S5 and visualized in Figure S8.The total BA for each dataset during the analysis period is given in Figure 7a, with the related box and whisker plots in Figure 7b.Despite the total area covering only 35% of the central/northern Siberia region and 30% of the North America region, southern Siberia sees a mean BA of 11.2 ± 4.46 (sd) Mha/year from 2001 to 2016.The years 2003 and 2008 were the largest fire years with 18.88 ± 10.61 and 22.88 ± 10.28 Mha BA on average across the datasets, respectively.However, unlike in North America and central/northern Siberia, southern Siberia has seen a decreasing tendency in the amount of BA per year, falling from 11.9 ± 5.56 Mha/year in 2001-2008 to 8.5 ± 3.25 Mha/year in the second half of the analysis period.The seasonality of the Talucci et al. polygon product can be compared to that of the MTBS/NBAC polygon product in North America, as it shows a tendency towards earlier fire occurrence as it does not take into account the spread of the fire between pixels, rather giving the date of the initial fire for the whole polygon.This can be somewhat improved by using the AvgDay value, which accounts for an increased proportion of late summer BA.However, both products show a single peak at approximately 35.4% in July following a constantly increasing number of fires over the previous four months and preceding a rapid two-month drop.

Southern Siberia
For southern Siberia the analysis is limited to the six global BA datasets (FireCCI51, MCD64A1, GABAM, GFED4s, GFED5, FireCCILT11) and MODIS active fires for the years 2001-2016.The results for the whole 2001-2020 period are shown in Table S5 and visualized in Figure S8.The total BA for each dataset during the analysis period is given in Figure 7a, with the related box and whisker plots in Figure 7b.Despite the total area covering only 35% of the central/northern Siberia region and 30% of the North America region, southern Siberia sees a mean BA of 11.2 ± 4.46 (sd) Mha/year from 2001 to 2016.The years 2003 and 2008 were the largest fire years with 18.88 ± 10.61 and 22.88 ± 10.28 Mha BA on average across the datasets, respectively.However, unlike in North America and central/northern Siberia, southern Siberia has seen a decreasing tendency in the amount of BA per year, falling from 11.9 ± 5.56 Mha/year in 2001-2008 to 8.5 ± 3.25 Mha/year in the second half of the analysis period.
Furthermore, there are large disparities between the datasets, even greater than in central/northern Siberia, with GFED5 showing over four times more BA than GABAM and GFED4s (Figure S9).These differences are further exacerbated during the smaller fire years; for instance, in 2002, 2004, and 2012, the BA in GFED5 is 530%, 758%, and 837% higher than GABAM, and 558%, 452%, and 518% higher than GFED4s, respectively.However, the largest difference between any dataset is in 2001, where GFED5 has over 10 times more BA than MCD64A1 (28.5 Mha vs. 2.82 Mha).
In southern Siberia the datasets can be split into three groups based on the statistical significance (Figure 7b and Table S6), as described in Appendix B. This time GFED4s produces the least amount of BA, 5.32 ± 2.94 Mha/year, whereas GFED5, with an average annual BA of 22.71 ± 7.91 Mha/year, produces more BA than all the others at a statistically significant level.Furthermore, there are large disparities between the datasets, even greater than in central/northern Siberia, with GFED5 showing over four times more BA than GABAM and GFED4s (Figure S9).These differences are further exacerbated during the smaller fire years; for instance, in 2002, 2004, and 2012, the BA in GFED5 is 530%, 758%, and 837% higher than GABAM, and 558%, 452%, and 518% higher than GFED4s, respectively.However, the largest difference between any dataset is in 2001, where GFED5 has over 10 times more BA than MCD64A1 (28.5 Mha vs. 2.82 Mha).
In southern Siberia the datasets can be split into three groups based on the statistical significance (Figure 7b and Table S6), as described in Appendix B. This time GFED4s produces the least amount of BA, 5.32 ± 2.94 Mha/year, whereas GFED5, with an average annual BA of 22.71 ± 7.91 Mha/year, produces more BA than all the others at a statistically significant level.Four of the six datasets where the Julian Day was provided (FireCCI51, GFED4s, GFED5, and MODIS active fires) agree on the seasonality of fires in southern Siberia (Figure 8).There is a sharp increase in the amount of BA from March to April, with a similar amount of BA produced in April and May (ca.35% each), before a sharp decrease to June.Unlike in North America and central/northern Siberia, a far smaller proportion of BA occurs during the summer months, with ca.6.2% happening during June, July, and August combined.There is then a second, smaller peak (ca.11.7%) in October.As in the other study regions, FireCCI51 shows a tendency towards later-occurring BA; however, this is not as notable here.The dataset shows an increase from April to May and has the lowest amount of BA produced between June and September, before the largest peak in October (14.1%).On the contrary, GFED5 shows that BA happens earlier in the season, showing a drop from 38.9% in April to 33.5% in May and the smallest peak in October (10.3%).
combined.There is then a second, smaller peak (ca.11.7%) in October.As in the other study regions, FireCCI51 shows a tendency towards later-occurring BA; however, this is not as notable here.The dataset shows an increase from April to May and has the lowest amount of BA produced between June and September, before the largest peak in October (14.1%).On the contrary, GFED5 shows that BA happens earlier in the season, showing a drop from 38.9% in April to 33.5% in May and the smallest peak in October (10.3%).MCD64A1 follows a more pronounced early-season BA occurrence tendency than the aforementioned datasets.Rather than a sudden increase/decrease between March and April/May and June, respectively, there is a single peak in April, with a greater proportion of BA produced in March (12.6%) and much less in May (21.2%).The summer and autumn BA occurrence patterns are similar to those of the other products although, as in spring, autumn BA in MCD64A1 happens slightly earlier, with September having ca.2% more BA; however, this difference is relatively minor.Conversely, the seasonality of FireCCILT11 shows a heavily weighted tendency towards later spring BA.Here just 23.8% of pixels are burned in April, with a substantial 57.0% of pixels burning in May, before a large drop to 2.8% in June and following similar tendencies to the other products to the end of the year.

Seasonal Patterns
Both North America and central/northern Siberia have broadly similar seasonal fire patterns, despite being composed of different tree types with distinct fire regimes.Both regions have their largest peak in July.However, while this is the only peak in North America, there is a second smaller, early-season peak in May in central/northern Siberia.This is likely due to the imperfect delineation between the two Siberia study regions, so the southern parts of central/northern Siberia are catching some of the spring fire regime from southern Siberia.The seasonality of fires in southern Siberia differs significantly from those in the more northerly study regions.The majority (ca.70%) of the BA here happens during April and May, with relatively few fires occurring during the summer months, before a second smaller peak (ca.11%) in October.Most precipitation in southern MCD64A1 follows a more pronounced early-season BA occurrence tendency than the aforementioned datasets.Rather than a sudden increase/decrease between March and April/May and June, respectively, there is a single peak in April, with a greater proportion of BA produced in March (12.6%) and much less in May (21.2%).The summer and autumn BA occurrence patterns are similar to those of the other products although, as in spring, autumn BA in MCD64A1 happens slightly earlier, with September having ca.2% more BA; however, this difference is relatively minor.Conversely, the seasonality of FireCCILT11 shows a heavily weighted tendency towards later spring BA.Here just 23.8% of pixels are burned in April, with a substantial 57.0% of pixels burning in May, before a large drop to 2.8% in June and following similar tendencies to the other products to the end of the year.

Seasonal Patterns
Both North America and central/northern Siberia have broadly similar seasonal fire patterns, despite being composed of different tree types with distinct fire regimes.Both regions have their largest peak in July.However, while this is the only peak in North America, there is a second smaller, early-season peak in May in central/northern Siberia.This is likely due to the imperfect delineation between the two Siberia study regions, so the southern parts of central/northern Siberia are catching some of the spring fire regime from southern Siberia.The seasonality of fires in southern Siberia differs significantly from those in the more northerly study regions.The majority (ca.70%) of the BA here happens during April and May, with relatively few fires occurring during the summer months, before a second smaller peak (ca.11%) in October.Most precipitation in southern Siberia occurs during July and August [29], which accounts for why there are fewer fires in the summer months here.The grass and shrubs, which grow quickly during this regenerative period, produce the ground-level fuel for repeat, low-moderate intensity fires in the 20-year period.

Advantages and Limitations of Different Methodologies
Compared to Siberia, the BA products show relatively good consistency in North America.This is partly since the official MTBS and NBAC databases are considered good representations of the total BA and can be used for validation, and partly because it is comparatively easier to detect BA from crown fires using remote sensing imagery due to the drastic change in pre-and post-fire canopy vegetation.In this regard, dNBR calculated from high resolution images from satellites such as Landsat can more precisely and successfully determine the extent of fires in North America, as with the MTBS and NBAC polygons and ABoVE-FED.
The products face significant challenges in determining the BA in both Siberian study regions, as the majority of fires are surface fires of low to moderate intensity [43,65].The canopy can obscure these fires from satellite imagery, and the lower energy output means that they are more difficult to detect using thermal sensors [42].In southern Siberia, the vast majority of the yearly fires occur during spring; these fires are often smaller in size but greater in quantity [92].The steppe then regenerates quickly during the rainy summer months with grass and shrubs growing between the burnt tree stands [47].As such, it is substantially more difficult for products using image composites at a low temporal resolution, such as annual, to determine where the fires took place, and this leads to large inconsistencies and variations between the datasets.The clearest example of these differences occurs in 2012: GFED4s has 0.77 Mha BA in central/northern Siberia and 3.49 Mha in southern Siberia, whereas GFED5 has 8.94 Mha and 21.58 Mha, respectively.Short-term composites of dNBR, such as with ABBA, again produce the most consistently strong results in Siberia; however, without reliable official data we cannot determine with any certainty the most accurate product.
Products with a greater omission error than commission error, such as FireCCI51, MCD64A1, and GABAM, will underestimate the total BA, and this has been shown in our analysis.For pixels to be declared burnt using the GABAM algorithm, the stringent threshold tests, using the surface reflectance of six bands alongside eight spectral indices, mean the total BA is likely to be a conservative estimate.ABoVE-FED was designed to roughly split omission and commission equally, so the amount of overestimated and underestimated BA is designed to be similar, thus producing a balanced total.
When considering which BA product to use, we suggest those with high spatial (≤250 m) and temporal (monthly) resolutions, and those which are made using short-term dNBR composites and include threshold tests.If applicable, products with roughly equal commission and omission errors and those which are designed for one particular region are also preferable.If there is concern over using a less reliable product, consider using another as validation.

Advantages and Limitations of Each BA Product
Even though FireCCI51 was masked to only show fires of 80% and greater confidence, it still performs well in the boreal regions, produces significantly more BA than MCD64A1 in southern Siberia and is more reliable than FireCCILT11.Alongside its relatively high spatial resolution, FireCCI51 is continually updated to the present month.Its availability as both a global pixel and gridded product allows for greater flexibility than other, single-type products.FireCCI51 is used in training both GABAM and FireCCILT11, demonstrating that it is considered versatile and reliable.MCD64A1 is also deemed a reliable product as it is used as the baseline dataset for GFED4s (Collection 5) and GFED5 (Collection 6).However, both versions have known limitations in boreal regions: C5 MCD64A1 underestimates in boreal Russia as fires detected by the MODIS hotspot algorithm are often too small to class the entire pixel as burned [124].Due to the masking of small lakes at high latitudes in North America, the BA in this region is significantly underestimated in C6 [71,125], which is why GFED5 has less BA in North America than GFED4s.
Despite the higher spatial resolution, the errors due to Landsat-7 scan line failings in GABAM for the years 2001-2013, which are not accounted for in the dataset methodology [72], are significant and users should be aware of its limitations in the study regions considered here for these years.The years 2014-2020 can be considered more reliable due to the introduction of Landsat-8 data from 2014 onwards, which supersedes data from Landsat-7.
It also appears that FireCCILT11 has significant limitations in representing Arcticboreal fire regimes, with the relatively low confidence level (30%) and spatial resolution (ca.5.5 km) negatively affecting the output.This is most evident with the fire seasonality, where FireCCILT11 significantly differs from all the other datasets in North America and southern Siberia.FireCCILT11 is also the only product which uses AVHRR data.AVHRR sensors have been in operation since 1981, and as such are fairly limited and have been known to degrade over time, producing more unreliable imagery with reduced quality compared to images from more modern satellites [126][127][128].
GFED5 shows a considerable improvement over GFED4s, and it performs well in both North America and central/northern Siberia, making it a reliable gridded product.However, when a dataset produces more BA than all the others at a statistically significant level, as with GFED5 in southern Siberia (Table S6), caution should be taken.
Region-specific products are, in theory, optimized for the local fire and vegetation characteristics.In North America, ABoVE-FED balances the omission and commission errors, and, by using the burnable fraction of a 500 m pixel derived from 30 m Landsat images, produces more accurate results than the comparable MCD64A1.ABBA performs similarly to ABoVE-FED, although ABoVE-FED provides data on the estimated day of burn for each pixel, whereas ABBA does not.The BA products in North America are compared to the official statistics from the 30 m MTBS and NBAC products, which are considered as reference BA.
In central/northern Siberia, however, we do not have a similar baseline from which we can compare the performance of all the products because of the absence of reliable official statistics.The Talucci et al. (2022) [63] dataset is an excellent representation of fire perimeters in central/northern Siberia at a high spatial resolution.However, the product does not always take into account unburned patches within the polygons, so it may overestimate the total BA.In this region, ABBA detects substantially more BA than most products, including the Russia-specific SRBA product, but less than GFED5, which, like Talucci et al. (2022) [63], may overestimate BA.SRBA produces less BA than both the global FireCCI51 and MCD64A1 datasets in central/northern Siberia, suggesting that the total BA is likely to be underestimated.
Active fire data, from MODIS or otherwise, are used in different ways to create BA products, but should not be taken as a definitive account of BA as they tend to overestimate small fires, but underestimate larger fires [124,129].In our study, with fires of low (less than 30%) confidence masked out, this is most notable in Alaska where the MODIS active fires vastly underestimate the total BA, especially in the large fire year of 2004, as also shown in Loboda et al. (2011) [130].The MODIS hotspot data are best used to determine which BA products perform well when considering the seasonality, as the first thermal anomaly for each pixel is determined twice per day (once for each pass of the Terra and Aqua satellites), thus modeling the spread of fire well.Active fire products from VIIRS have been shown to perform better than those from MODIS [131,132]; however, we did not include these due to the relatively short time series available.

Study Limitations
This study took measures where necessary, as described in Section 2.3.1, to ensure that any re-occurrences of fire events within a single year were dealt with.Another limitation is that the Julian Day values for the polygon products are not directly comparable to the Julian Day values for the pixel products.The pixel products show the dynamic nature of each fire as different values are given for each pixel within the same fire event, based upon how the fire spreads.For the polygon products the Julian Day value is assigned for the entire polygon, meaning that the seasonality of the fires will be biased towards early season fires, as shown especially with the Talucci et al. fire perimeter product.
As shown in Table 1, MCD64A1 is analyzed with other products across several other studies, with either active fire data or high-resolution imagery used as validation.Zhang et al. (2018) [61] look at how GFED4, GFED4s, and MCD64A1 (C6) compare with both MODIS active fire data and higher-resolution imagery from VIIRS, Landsat, and Sentinel-2 in East China (2015-2016) and North-West India (2016 only) to determine which dataset best represents agricultural burning in these regions.Fires in China, specifically the mountainous area of Northwest Yunnan, are also the focus in Fornacca et al. (2017) [56], which looks at MCD45A1 (C5.1),MCD64A1 (C6), and FireCCI41 alongside MODIS active fire data and Landsat reference imagery for two sample years (2006 and 2009) to determine which dataset performs best in representing small fires in this area of steep terrain.MCD64A1 (C6), GFED4, and GFED4s are also compared in Mediterranean Europe, this time with FireCCI51 and local data from the European Forest Fire Information System for the years 2001-2011 [59].Forty-four areas across Northern Hemisphere Africa and South America were used to validate the output from MCD45A1 (C5), MCD64A1 (C5 and C6), and FireCCI (4.1 and 5.0) for the years 2007-2008 in Valencia et al. (2020) [60].Staying in South America, parts of the Brazilian Amazon biome have also been studied in Pessôa et al. (2020) [58] and Dutra et al. (2023) [55].The former considers the global MCD64A1 (C6), FireCCI50, and GABAM products against region-specific data from the Tropical Ecosystems and Environmental Sciences lab for fires during the summer and autumn months of 2015, whereas the latter compares MCD64A1 (C6) and GABAM against the Global Wildfire Information System [133] and the regional MapBiomas Fogo product [134] for fires in 2019.Humber et al. (2019) [57] is the only study to use a similar methodology to the present work.However, the paper considers the total BA globally, split into several ecoregions as defined by Giglio et al. (2006) [108] and van der Werf et al. (2006) [135], for the years 2005-2011.MCD45A1 (C5), MCD64A1 (C6), and FireCCI41 are compared with MODIS active fires and also the Copernicus Burnt Area product.It is difficult to directly compare our results to those from previous works due to the use of different versions of datasets (e.g., FireCCI51 in this work vs. FireCCI41 in Humber et al. (2019) [57]) and the products perform differently in boreal forests compared to other regions.

Conclusions
To conclude, we conducted what we believe to be the first comprehensive analysis of 12 BA products alongside MODIS active fire data across North America, central/northern Siberia, and southern Siberia.We found that there is relatively little variability between the datasets in North America both in terms of total BA and seasonality, whereas in both central/northern and southern Siberia there are large disparities between the products.These disparities occur due to the different methodologies used to determine BA, with dNBR of short-term composites from Landsat images used alongside hotspot data seemingly the most consistently successful.Differences between pixel, gridded, and polygon products should be considered when choosing a BA dataset.We suggest that FireCCI51 and MCD64A1 are relatively reliable and consistent global pixel products, and GFED5 can be considered a reliable global gridded product.On the other hand, we stress caution using GABAM, especially for the years 2001-2013 in these regions.We also highlight ABoVE-FED and ABBA for use specifically in North America, and Talucci et al. in central/northern Siberia, due to their detection of smaller fires which are often missed by global products.Future work could consider extending the analysis into the 2020s and including newly-developed datasets, such as FireCCIS311 [136], which uses Sentinel-3 images alongside VIIRS active fire data.Furthermore, projections of future fire occurrence could be made under different climate change scenarios to determine the potential risks of fire to indigenous communities and the impact it could have in the wider Arctic-boreal zone.those with p-values less than 0.05.The outputs from this test for each region are given in Tables S2, S4 and S6.
To conduct these tests by hand, the test statistic t is calculated using Equation (A2), and the degrees of freedom value df is calculated using Equation (A3).where x is the mean, s is the standard deviation, and n is the sample size for each pair of datasets i, j to be compared.df is rounded down to the nearest integer and the p-value is then found using statistical tables.

36 Figure 1 .
Figure 1.Map of study regions-North America (red), central/northern Siberia (blue), and southern Siberia (purple).The internal borders of each ecoregion in North America and central/northern Siberia are visible.

Figure 1 .
Figure 1.Map of study regions-North America (red), central/northern Siberia (blue), and southern Siberia (purple).The internal borders of each ecoregion in North America and central/northern Siberia are visible.

Figure 2 .
Figure 2. Total BA (a) and box and whisker plot (b) for each dataset for the analysis period (2002-2016) in North America.The whiskers extend to 1.5 * Inter-Quartile Range.The mean of each dataset is given as the green triangle.

Figure 2 .
Figure 2. Total BA (a) and box and whisker plot (b) for each dataset for the analysis period (2002-2016) in North America.The whiskers extend to 1.5 * Inter-Quartile Range.The mean of each dataset is given as the green triangle.

Figure 4 .
Figure 4. Total BA (a) and box and whisker plot (b) for each dataset for the analysis period (2006-2016) in central/northern Siberia.The whiskers extend to 1.5 * Inter-Quartile Range.The mean of each dataset is given as the green triangle."Talucci et al." refers to the Talucci et al. (2022) fire perimeter product [63].

Figure 4 .
Figure 4. Total BA (a) and box and whisker plot (b) for each dataset for the analysis period (2006-2016) in central/northern Siberia.The whiskers extend to 1.5 * Inter-Quartile Range.The mean of each dataset is given as the green triangle."Talucci et al." refers to the Talucci et al. (2022) fire perimeter product [63].

Figure 5 .
Figure 5. (a) The "stripe" and "box" effects in GABAM shown in 2001 and (b) scan lines appearing as missing patches through burned areas in 2003 in the central/northern Siberia region.When considering the seasonal patterns over the 2001-2016 period, as shown in Figure6, the datasets show further disagreement in central/northern Siberia compared to North America.GFED4s, GFED5, and MODIS active fires show a constant increase to an initial peak (ca.16%) in May, followed by a drop to 12.5% in June, a second, larger peak ranging from 28 to 38% in July, and a steady decrease through August to ca. 5.5% in September, before dropping below 1.5% in October.

Figure 5 .
Figure 5. (a) The "stripe" and "box" effects in GABAM shown in 2001 and (b) scan lines appearing as missing patches through burned areas in 2003 in the central/northern Siberia region.

Figure 7 .
Figure 7.Total BA (a) and box and whisker plot (b) for each dataset for the analysis period (2001-2016) in southern Siberia.The whiskers extend to 1.5 * Inter-Quartile Range.The mean of each dataset is given as the green triangle.

Figure 7 .
Figure 7.Total BA (a) and box and whisker plot (b) for each dataset for the analysis period (2001-2016) in southern Siberia.The whiskers extend to 1.5 * Inter-Quartile Range.The mean of each dataset is given as the green triangle.

Table 1 .
A list of previous relevant literature with associated study areas, products used, plus years of interest.

Table 2 .
Information about the BA datasets used in this study.