Future projections of Siberian wildfire and aerosol emissions

Wildfires are among the most influential disturbances affecting ecosystem structure and biogeochemical cycles in Siberia. Therefore, accurate fire modeling via dynamic global vegetation models is important for predicting greenhouse gas emissions and other biomass-burning emissions to understand changes in biogeochemical cycles. We integrated the widely used SPread and InTensity of FIRE (SPITFIRE) fire module into the spatially explicit individual-based dynamic global vegetation model (SEIB-DGVM) to improve the accuracy of fire predictions and then simulated future fire regimes to better understand their impacts. The model can reproduce the spatiotemporal variation in biomass, fire intensity, and fire-related emissions well compared to the recent satellite-based estimations: aboveground biomass (R²=0.847, RMSE =18.3 Mg ha⁻¹), burned fraction (R²=0.75, RMSE=0.01), burned area (R²=0.609, RMSE =690 ha), dry-matter emissions (R²=0.624, RMSE =0.01 kg DM m⁻²; dry matter), and CO₂ emissions (R²=0.705, RMSE =6.79 Tg). We then predicted that all of the 33 fire-related gas and aerosol emissions would increase in the future due to the enhanced amount of litter as fuel load from increasing forest biomass production under climate forcing of four Representative Concentration Pathways: RCP8.5, RCP6.0, RCP4.5, and RCP2.6. The simulation under RCP8.5 showed that the CO₂, CO, PM_2.5, total particulate matter (TPM), and total particulate carbon (TPC) emissions in Siberia in the present period (2000–2020) will increase relatively by 189.66±6.55, 15.18±0.52, 2.47±0.09, 1.87±0.06, and 1.30±0.04 Tg species yr⁻¹, respectively, in the future period (2081–2100) and the number of burned trees will increase by 100 %, resulting in a 385.19±40.4 g C m⁻² yr⁻¹ loss of net primary production (NPP). Another key finding is that the higher litter moisture by higher precipitation would relatively suppress the increment of fire-related emissions; thus the simulation under RCP8.5 showed the lowest emissions among RCPs. Our study offers insights into future fire regimes and development strategies for enhancing regional resilience and for mitigating the broader environmental consequences of fire activity in Siberia.

1 Introduction

Fires are among the most significant disturbances affecting biogeochemical cycles, atmospheric chemistry, the carbon cycle, and ecosystem structure and function worldwide (Pickett et al., 1999). Wildfires are also the dominant climate-driven disturbance agent in boreal forests (Goldammer and Furyaev, 1996; Shorohova et al., 2011; De Groot et al., 2013), shaping major forest cover in Russia (Krylov et al., 2014) and rapidly increasing the burned area and emission intensity in Canada and Alaska (Zheng et al., 2021). Fires influence vegetation dynamics by allowing plants to adapt to fire regimes, influencing vegetation productivity, litter, and fuel load (Cochrane, 2003; Bergeron et al., 2004; Whelan, 2009). The intensity and frequency of large-scale boreal forest fires are expected to increase in the future due to increased global temperatures, drier conditions, and longer fire seasons, which will cause more emissions from biomass burning (Flannigan et al., 2009; Gauthier et al., 2015) and human activity using fire for land management (e.g., use of fire as a tool in the deforestation process) (Hantson et al., 2016; Archibald et al., 2013; Morton et al., 2008). Globally, from 2000 to 2019, satellites detected a decrease in the burned area of grassland, while there was a slight increase in the area of forest fires in Russia (Zheng et al., 2021). Central and southern European countries, such as France, Spain, Portugal, and Greece, are already experiencing larger and more devastating fires (Carnicer et al., 2022). Not only large fires but also small fires have a significant impact: areas burned by small fires contributed 35 % to the total burned area, from 345 to 464 Mha yr⁻¹, and related carbon emissions increased from 1.9 to 2.5 Pg C yr⁻¹ from 2001–2010 (Randerson et al., 2012). This finding is in line with current studies reporting that the global mean CO₂ emission intensity has increased by 0.9±0.9 % yr⁻¹ from 2000 to 2019 (Zheng et al., 2021) and that the fire weather index (FWI) reached levels above 30, corresponding to high, very high, and extreme levels of fire frequency, causing CO₂ emissions to increase in Europe since 1980 (Carnicer et al., 2022).

Forest fires are important ecological factors that influence both the establishment and succession of vegetation (Abaimov and Sofronov, 1996). Climate-driven large fires are responsible for rapid changes in vegetation (Cleve and Viereck, 1981), soil properties (Pastor and Post, 1986; Pellegrini et al., 2021), biogeochemical cycling, microclimate, forest ecosystems (Crutzen and Goldammer, 1993), productivity, stability, and many other ecological properties (Melillo et al., 1993). Forest fires also indirectly affect vegetation dynamics by increasing CO₂ levels in the atmosphere (Seiler and Crutzen, 1980; Nguyen and Wooster, 2020), as CO₂ is one of the primary products of biomass combustion and is emitted in all phases of fire (ignition, flaming, glowing, pyrolysis, and extinction) (Andreae and Merlet, 2001), with the flaming phase leading to emissions (Lobert et al., 1991; Ward and Hardy, 1991). Thus, it is challenging to estimate CO₂ emissions because they are generated in large quantities during biomass combustion and because of the different emission timelines produced during each combustion stage. Increasing atmospheric CO₂ concentrations alter the global carbon cycle by causing global warming (Van Der Werf et al., 2006, 2010, 2017; Neto et al., 2009; Kaiser et al., 2012; Lin et al., 2013), and the resulting global warming is expected to intensify extreme fire seasons, leading to further surges in carbon emissions that significantly contribute to the global burden of greenhouse gases (fire–climate feedbacks) (Bowman et al., 2009). Therefore, accurate modeling of future wildfires and their emissions is required to understand the associated risks.

Boreal vegetation stores 17 % of the world's carbon yet encompasses almost 30 % of all terrestrial carbon stocks (Kasischke, 2000; Gauthier et al., 2015), with two-thirds being located in Siberia, Russia (Shvidenko and Nilsson, 2003). In Siberia, burned-biomass emissions approached 0.4 Gt C yr⁻¹ in 2021, 3 times the average value between 1997 and 2020, according to GFED4s (Friedlingstein et al., 2020). Kharuk et al. (2022) also stated that the decadal frequency of wildfires tripled between 2001–2010 and 2011–2020. Catastrophic boreal forest fires are expected to continue to increase in the future due to increased global temperature, drier conditions, and longer fire seasons, and these fires will increase the severity and emissions produced from biomass burning (Flannigan et al., 2009). Burning vegetation is a major source of black carbon (BC), carbon monoxide (CO) (Forster et al., 2018), and particulate matter (PM) (Reddington et al., 2016). According to records from the Copernicus Atmosphere Monitoring Service (CAMS), Russia experienced a drastic increase to 8 Mt (metric megatons) in PM_2.5 emissions in 2021, which is 78 % higher than the average level between 2004 and 2021 (4.5 Mt) (Romanov et al., 2022). Furthermore, an increase in atmospheric emissions negatively affects the climate by contributing to global warming and climate change (Randerson et al., 2006; Westerling et al., 2006; Bowman et al., 2009) and affects weather systems by modulating solar radiation and cloud properties (Schultz et al., 2008).

Understanding how long-term climate change, fire regimes, and forest vegetation interact under multiple climate scenarios is critical for forecasting forest succession trends (Clark and Richard, 1996). Modeling of fire regimes using dynamic global vegetation models (DGVMs) is a key approach to analyzing these factors. However, including interactive fire disturbances in vegetation models is critical for accurately simulating vegetation dynamics (Thonicke et al., 2001). A well-structured process-based fire module can accurately assess fire activity, consumed biomass due to fire, and biomass-burning emissions. The assessment of each fire-related variable is interconnected with another variable, so the module must be well constructed because the amount of consumed biomass during forest fires can vary significantly. Several factors affect burned biomass, such as spatial and temporal variations in the burned area based on the ignition factors, quantity and quality of the fuels available, and vegetation or plant functional type (PFT); additionally, every PFT reacts differently to fire disturbance (Cramer et al., 2001; Ito, 2011). Since the first global fire models were integrated into dynamic global vegetation models (DGVMs) 2 decades ago, the variety and complexity of fire models have expanded (Hantson et al., 2016). The Fire Modeling Intercomparison Project (FireMIP) compared 11 current fire models by structure and simulation protocols, using a benchmarking system to evaluate the models (Rabin et al., 2017). The results indicate that models that explicitly distinguish ignition factors, such as lightning and human-caused “ignition events”, as well as physical properties and processes that determine fire spread and intensity by plant functional type (PFT), performed better. One such fire module is SPITFIRE (SPread and InTensity of FIRE; an upgrade of GlobFIRM) (Thonicke et al., 2010), which has been used in several DGVMs: LPJ-GUESS–SPITFIRE, ORCHIDEE–SPITFIRE, JSBACH–SPITFIRE, and LPJ-LMfire. In this study, we integrated the SPITFIRE fire module into the spatially explicit individual-based dynamic global vegetation model (SEIB-DGVM) to predict fire, vegetation, and burned-biomass emission variables in Siberia in the future. We selected SEIB-DGVM because of its high-quality biogeochemical model coupled with a three-dimensional representation of forest structure where individual trees compete for light and space (Sato et al., 2007). The SEIB-DGVM processes physical, physiological, and vegetation dynamics and was previously used for reconstructing the geographical distributions of fundamental plant productivity properties (Sato et al., 2020), evaluating the geographic and environmental heterogeneity of larch forests with a special focus on topography (Sato and Kobayashi, 2018), and assessing the impacts of global warming on Siberian larch forests and their interactions with vegetation dynamics and thermohydrology (Sato et al., 2016). SEIB-DGVM accurately simulates forest ecology after typhoon disturbances (Wu et al., 2019), nonstructural carbohydrate dynamics (Ninomiya et al., 2023), and masting in a temperate forest (Végh and Kato, 2024).

The original fire module of SEIB-DGVM is Glob-FIRM (Thonicke et al., 2001), which has several limitations; for example, human-changed fire regimes and other land use impacts are not considered (Thonicke et al., 2001). In addition, GlobFIRM derives the burned fractional area of a grid cell from the simulated length of the fire season and from the minimum annual fuel load; this method does not specify ignition sources and assumes a constant fire-induced mortality rate for each plant functional type (PFT) (Thonicke et al., 2010). To improve the fire simulations with SEIB-DGVM, we replaced its fire module with the SPITFIRE model (Thonicke et al., 2010) by adding complete ignition equations (human and lightning effects, etc.). The module included a calculation mechanism for trace gas and aerosol emissions (Andreae and Merlet, 2001) and was adjusted to produce monthly outputs for all variables in SEIB-DGVM. These improvements allowed us to simulate fire activity and aboveground biomass dynamics and spatiotemporally assess the projected burned biomass and its emissions for the 21st century in Siberia under Representative Concentration Pathways (RCPs).

2 Methods

2.1 Study sites

Boreal forests represent the largest forest biome and one-third of global forest cover (De Groot et al., 2013) and play an important role in the atmosphere–land interactions of the global climate system (Randerson et al., 2006; Bonan, 2008). Geographically, boreal forests are found in Canada, Alaska, and Siberia, of which Siberia has the largest forested area. Siberia is largely covered by deciduous needleleaf conifers (Fig. 1), which consist mostly of the larch species Larix sibirica, L. decidua, and L. dahurica (Abaimov et al., 1998), which are categorized as pyrophytic species, meaning that they require periodic fires to persist on the landscape (Kharuk et al., 2011). The Siberian land cover has changed very little over the last century (Ivanov et al., 2022), and the boreal forest covers approximately $>\mathrm{15}\times {\mathrm{10}}^{\mathrm{6}}$ km², containing a large amount of carbon that is comparable to the combined carbon storage in tropical and temperate forests (Dixon et al., 1994; Kasischke, 2000).

The main external factors affecting Siberian boreal forests are fires and climate change (Goldammer and Furyaev, 1996; Shorohova et al., 2009). Climate change increased the frequency of forest fires, which in turn amplified the impacts of climate change locally. In the Arctic, a rapid warming trend has been observed, and the increase in temperature over the last 20 years of the 20th century was 2 to 3 times higher than the global average, while in the first 20 years of the 21st century, it exceeded 4 times (Chylek et al., 2022). This enormous increase in temperature in Siberia affected the duration and speed of snowmelt and accelerated the thawing of carbon-rich permafrost (Natali et al., 2019; Schuur et al., 2015; Nitzbon et al., 2020), which results in drier ground cover, an increased frequency of wildfires, longer fire seasons, and increased ignition sources (Kharuk et al., 2022). These changes may result in a new climate state in which heat waves as well as the associated the occurrence of wildfires may become routine and more severe (Hantemirov et al., 2022; Landrum and Holland, 2020). Produced emissions from thawing permafrost and from wildfires are likely to feed into the global carbon cycle's feedback on climate change (Schuur et al., 2015), triggering further warming trends globally (Schimel et al., 2001; Kharuk et al., 2011; Krylov et al., 2014).

Figure 1 Study site (black rectangle: 60–180° E, 45–80° N). Green and brown indicate the forest types in Siberia as provided by the Global Land Cover dataset (GLC 2000) for northern Eurasia (Bartalev et al., 2003). Grey indicates other vegetation types in the Siberian area as provided by the Database of Global Administrative Areas (GADM).

2.2 Improved fire module principles

We improved the SEIB-DGVM fire module by replacing the Glob-FIRM (Thonicke et al., 2001) with the SPITFIRE model (Thonicke et al., 2010). First, we added two new input variables for fire ignition: population and lightning data. Second, we incorporated the complete SPITFIRE equation (Thonicke et al., 2010), which included new variables, PFT parameters, and local parameters, and improved the output to be able to be produced on a monthly scale (Fig. 2). The variable integration between the default and improved fire models requires several parameter-specific PFTs (Table 1).

The default SEIB-DGVM model uses annual time steps for vegetation dynamics and disturbance, which we improved to monthly time step outputs. The fraction of individual trees killed by a fire depends on PFT fire resistance (M₃, Table 1). All grass leaf biomass, all dead and living tree leaf biomass, half of the dead tree trunk biomass, and half of the litter pool are released into the atmosphere as CO₂ during a fire, while the dead tree's residual biomass is converted into litter. In reaction to fire, all deciduous PFTs convert their phenology phase to dormancy, and if the stock resource of grass PFTs (gmass stock) does not meet the minimal value (50 g DM m⁻²; dry matter) following a fire, the deficit is supplemented from litter (Sato et al., 2007). Furthermore, related to the fire–vegetation relationships, for herbaceous PFTs, both belowground and storage biomass are preserved after a wildfire and used for the recovery of aboveground biomass. During this recovery period, herbaceous PFTs work on producing aboveground biomass while reducing their storage biomass, thus increasing the allocation ratio to aboveground biomass in the post-fire phase. For woody PFTs, fire only gives the option for individual trees to either die or survive. The surviving trees only lose their foliage biomass. As the foliage is lost, fine-root biomass becomes unnecessary, leading to its rapid loss due to its fast turnover rate. In the spring following a fire, surviving trees convert storage resources into foliage and fine-root biomass. The new net primary production (NPP) from the newly formed foliage first prioritizes the recovery of leaves and fine roots. Therefore, fires increase the allocation ratio to the foliage and fine roots in surviving woody plants.

The basic equation of fire disturbance is the area burned, which we adjusted with the SPITFIRE equation (Thonicke et al., 2010) by including the fire probability and area of the grid cell:

$$\begin{array}{}\text{(1)}& A_{\mathrm{b}}=P_{\mathrm{b}}\times A,\end{array}$$

where A_b is the area burned in a grid cell per month (ha per month), P_b is the product of the probability of fire per month at any point inside the grid cell (per month), and A is the area of the grid cell (ha). P_b is the fire probability and is the product of the fuel load (litter + aboveground biomass) and its moisture factor. We used the same P_b mechanism as that of the default fire module, where if the fuel load satisfies the minimum fuel threshold (200 g C m⁻²), random fires can occur at any point location inside the grid cells. In this improvement, P_b was modified by considering the ignition event E(n_ig) (ha⁻¹) alongside anthropogenic- (human population density) and natural-ignition (lightning strikes) possibilities, the fire danger index (FDI), and the mean fire area $\stackrel{\mathrm{‾}}{a}f$(ha). Thus, Eq. (2) can be represented as follows:

$$\begin{array}{}\text{(2)}& A_{\mathrm{b}}=E\left(n_{\mathrm{ig}}\right)\times \mathrm{FDI}\times \stackrel{\mathrm{‾}}{a}f\times A.\end{array}$$

Technically, the SEIB-DGVM simulation of each grid cell is carried out independently among the surrounding grid cell, so the fire cannot spread to other grid cell without those grid cell meeting the ignition requirements (fuel load and fuel moisture).

After all variables in the SPITFIRE fire module were integrated, we added the trace gas and aerosol emission calculation process to the model. Trace gas and aerosol emission estimation are line with the Fire Modeling Intercomparison Project (FireMIP) protocols (Li et al., 2019), the comprehensive study comparison of nine dynamic global vegetation models (DGVMs), and produced an important estimation for long-term and large-scale fire emissions. Using the FireMIP protocol reference, SEIB-DGVM SPITFIRE improved to output PFT-level fire emissions.

Trace gas and aerosol emissions are the result of the total amount of burned biomass, the sum of dead and live fuel consumption as the result of surface fire, and crown scorch. Trace gas emissions are estimated based on fire carbon emissions, vegetation characteristics from DGVMs, and fire emission factors. Fire emissions of trace gas and aerosol for each species i and PFT j, E_i,j (g species m⁻²), are estimated based on Andreae and Merlet (2001):

$$\begin{array}{}\text{(3)}& E_{i,j}={\mathrm{EF}}_{i,j}\times {\displaystyle \frac{{\mathrm{CE}}_j}{C}},\end{array}$$

where EF_i,j is the PFT-specific emission factor (g species kg⁻¹ DM); CE_j is the combusted biomass of PFT_j due to the fire (g C m⁻²); and C is the unit conversion factor from carbon to dry matter, $C=\mathrm{0.5}\times {\mathrm{10}}^{\mathrm{3}}$ g C kg⁻¹ DM. The emission factors (EFs) used in this study are based on Andreae and Merlet (2001) and the updated pyrogenic emission species by various types of biomass burning (Andreae, 2019) (Table S2 in the Supplement).

DGVMs generally simulate vegetation as a combination of PFTs in a given grid location to represent plant function at a global scale, instead of land cover types (Li et al., 2019). In this, we classified the PFTs with the land cover types (LCTs) to integrate the emission factors of each LCT for trace gas and aerosol emission estimation processes. The BoNE (boreal needleleaf evergreen), BoNS (boreal needleleaf summergreen), and BoBS (boreal broadleaf summergreen) PFTs are classified as boreal forest LCTs. Other PFTs have been integrated with LCTs but are not listed in Table 1, as this study only covers boreal forest. The other integrated PFTs include TrBE (tropical broadleaf evergreen), TrBR (tropical broadleaf raingreen), TeNE (temperate needleleaf evergreen), TeBE (temperate broadleaf evergreen), TeBS (temperate broadleaf summergreen), TeH (temperate herbaceous C3 grass), and TrH (tropical herbaceous C4 grass). Related to the PFTs abbreviations, Tr represents tropical forest and Te represents temperate forest. SEIB-DGVM did not classify crop PFTs, so cropland LCTs will not be used.

Further changes to the input and output of the new SEIB-DGVM SPITFIRE are shown in Appendix A1, while Appendix A2 summarizes the improvement processes represented in this study, which can be classified into two groups: disturbance and biogeochemical dynamics. Appendix A3 lists the symbols used in the model's equations. Detailed information about the integration of the SPITFIRE module in SEIB-DGVM, which includes the improvement and adjustment of all the variables and the main important variables, such as the ignition events E(n_ig), fire danger index (FDI), mean fire area $\stackrel{\mathrm{‾}}{a}f$, fuel moisture content, rate of spread, fire fraction and intensity, fire damage to plants, and trace gas and aerosol emissions, is provided in the Supplement (Sect. S2.2.1–S2.2.7).

Figure 2 SEIB-DGVM SPITFIRE system diagram. In describing the improvements (SPITFIRE), the interaction with the previous fire module (Glob-FIRM) is shown. All original SPITFIRE variables were integrated: ignition factor (lightning and population), PFT parameters, and other fire-related variables. In addition to the default annual output, the improved module has monthly outputs of all variables depending on the user needs. For definitions of the abbreviations, refer to Appendix A3.

Table 1 SEIB-DGVM SPITFIRE PFT-specific (plant functional type) model parameter values and their attribution to LCTs. This table was modified from Thonicke et al. (2010).

PFTs attributed to land cover types (LCTs) are needed to classify the fire emission factor (EF) (Table S2) to estimate trace gas and aerosol emissions (Andreae and Merlet, 2001).

2.3 Model calibration

We calibrate the improved model using all of the benchmark datasets (Table 3). The calibration process is done sequentially for all of the major variables, burned fraction, burned area, dry matter, aboveground biomass, burned-biomass emissions, and the forest ecology variables (Fig. 3). This calibration process involves comparing the average value of the output variable with the corresponding variable from the benchmark dataset, ensuring that both are aligned in terms of spatiotemporal resolution. The process is sequential because one variable is used for the calculation of another variable (such as the burned fraction and burned area affecting aboveground biomass, forest structure, dry matter, and emissions). One calibration process is performed with multiple iterations until the output variable has numerical values and a spatial distribution similar to the benchmark data, and the process is repeated for other variables once the previous variable has been calibrated.

2.4 Model application

The original SEIB-DGVM utilizes three computational time steps: a daily time step for all physical and physiological processes except for soil decomposition and tree growth, a monthly time step for soil decomposition and tree growth, and an annual time step for vegetation dynamics and fire disturbance (Sato et al., 2007). In this study, we improved the fire module to calculate natural- and anthropogenic-fire-ignition factors (based on lightning flashes and population density) and adjusted it to produce monthly outputs using temporal-resolution statistical downscaling methods with user-defined weighted monthly parameters (Table 2). The annual average ignition factor variables (population density and lightning flash rate) were used consistently throughout all simulation phases.

We ran the improved model (SEIB-DGVM SPITFIRE) and the default model (SEIB-DGVM GlobFIRM) under the same protocols to equally compare and assess their fire products (Fig. S3 in the Supplement). Simulations were run in three phases (spin-up, historical, and future), and the simulation was run with the fire mode on and fire mode off to compare and assess the vegetation products during fire, and also each phase was replicated five times to minimize bias due to random variables in the tree morality¹ (see reference marks in Fig. 3). The model was run in three phases²: (1) a 1000-year spin-up phase to bring the soil and vegetation carbon pools into equilibrium with the climate using daily baseline CRU TS3.22 (Climatic Research Unit Time-Series) climate data; (2) a 156-year historical phase also using daily baseline CRU TS3.22 climate data and spin-up simulation results as inputs; and (3) a 95-year future phase using daily MIROC5 output submitted for the IPCC Fifth Assessment Report (AR5) (MirocAR5 base V3) with RCP8.5, RCP6.0, RCP4.5, and RCP2.6 climate data and historical simulation results as inputs (Fig. 3). The MirocAR5 base V3 dataset has been bias-corrected with CRU TS3.22 climate data, so using these two datasets consecutively in spin-up, historical, and future simulations ensures the harmony of the input climate data. Five different types of RCP scenario climate data were used to determine the impact of fire and climate on forest structure and their interactions.

In the previous SEIB-DGVM study, a 2000-year spin-up was needed to obtain the convergence amount of soil organic matter (Sato et al., 2010). However, we have conducted preliminary simulations with the same study area by setting the spin-up years to 1000 and 2000 years. We confirmed that the outputs of the 1000-year and 2000-year spin-up simulations were very similar; thus, the 1000-year spin-up was enough to reach carbon stock equilibrium. This parameter setting is also in line with the simulation settings in other SEIB-DGVM studies: Sato et al. (2007) performed a 1000-year spin-up and combined it with all of the simulation phases to extract general trends of post-fire succession. Another study by Arakida et al. (2021) also confirmed that a spin-up period of 100 years was sufficient for the equilibrium of the LAI (leaf area index), aboveground biomass, and GPP (gross primary production) at all the study sites in Siberia.

In addition, we have a verification stage³ to ensure that the new input data can be read, produced, and processed properly (Rabin et al., 2017). Then, we calibrate all of the major emissions individually and sequentially with the benchmark dataset because each variable affects other variables and we need to ensure the final output is comparable with the benchmark datasets⁴. After verifying that the new module was incorporated seamlessly, we validated the model outputs (fire, vegetation, and emission variables) using the GFED4, GFED4s, ESA Biomass Climate Change Initiative (CCI), and Global Biomass Burning Emissions Inventory (GBEI) benchmark datasets⁵.

Figure 3 Workflow of improving the SEIB-DGVM fire module.

Table 2 SEIB-DGVM SPITFIRE input data descriptions.

2.5 Model benchmarks

A common method for validating the outputs of dynamic global vegetation models (DGVMs) is to use satellite-based product datasets. For instance, direct observations of global fire occurrence by satellite-borne sensors can detect active fires, fire radiative power, and burned areas, and these observations have been available since the 1990s (Mouillot et al., 2014). The Fire Modeling Intercomparison Project (FireMIP) also used the satellite-based product database as a benchmark to evaluate the model simulation (Rabin et al., 2017; Li et al., 2019).

In the last few decades, several global biomass-burning emission datasets based on the detection of the burning area and fire radiative energy have been developed and used for many purposes, such as global climate and vegetation modeling, together with environmental, health, and security assessments (Ichoku et al., 2008; Mouillot et al., 2014). Although fire-related observation datasets are available and globally accessible, they have relatively large uncertainties and are poorly constrained, especially in models at the global and regional levels (Liousse et al., 2010; Petrenko et al., 2012, 2017; Bond et al., 2013; Zhang et al., 2014; Pan et al., 2015; Pereira et al., 2016).

Pan et al. (2020) reported that this uncertainty could be caused by various measurement and/or analysis processes, including the detection of fire or burned areas, retrieval of fire radiative power, emission factor information, biome type, burning stage, and fuel consumption estimation. The emission factor (EF) is considered an important factor for obtaining specific gaseous or particulate species of smoke emitted from burned dry matter in all major burned-biomass (BB) emission datasets. Some EFs originate from laboratory experiments where fuel samples are burned in combustion chambers (Christian et al., 2003; Freeborn et al., 2008), whereas others originate from large-scale, open biomass-burning and wildfire experiments. The combustion properties might differ greatly between these two categories; e.g., because of personnel security and other logistical considerations, some EF measurement locations are often not close enough to the biomass-burning source (Aurell et al., 2019). Another factor is the biome type, which affects the scaling factor of the emission coefficient for the FRP-based (fire radiative power) BB datasets (GFAS, Global Fire Assimilation System; FEER, Fire Energetics and Emissions Research; and QFED, Quick Fire Emissions Dataset). The emission factors of all BB datasets were assigned based on the type of biome, and most of the examined BB datasets had different definitions of major biome types, so uncertainty might be present at certain levels (Pan et al., 2020).

We validated the improved SEIB-DGVM fire module products using the burned-area (GFED4) and burned-fraction (GFED4s) datasets, corresponding to the model's output. These datasets have higher resolutions than other burned-area-based datasets, and all of the uncertainty probabilities regarding the selected database described by Pan et al. (2020) were adjusted with our model configurations. We used the emission factor (EF) from Andreae and Merlet (2001) with the latest update from Andreae (2019) and integrated the plant functional type (PFT) model with the land cover types (LCTs) used in the EF (Tables 1 and S2).

Furthermore, fire models should be evaluated together with their associated vegetation models because the former might produce burned areas perfectly but incorrectly simulate aboveground biomass (AGB) patterns. Fire products depend on AGB availability, and fire also affects AGB availability and succession after forest fires. Thus, to ensure that the model conducted correct assessments, we evaluated the aboveground biomass variable using the ESA Biomass Climate Change Initiative dataset (Table 3). The AGB data from the ESA Biomass Climate Change Initiative (CCI) v3 (2010, 2017, and 2018) include high-quality data with a large resolution of 100 m×100 m obtained from multiple remote sensing observations collected around the year 2010 (Santoro et al., 2021), making them suitable for validating our improved model product.

Overall, we validated the model spatially and numerically at the Siberian level and in smaller regions to determine the performance of the model in many points of view (spatial, numeric, wide, and small region). We classified Siberia into three regions: west region (60–90° E, 45–80° N), central region (90–120° E, 45–80° N), and east region (120–180° E, 45–80° N) (Fig. S12).

Table 3 Description of the observational datasets used for model evaluation.

3 Results

3.1 Improved model validation

3.1.1 Fire products

We compared the annual average distribution patterns of the burned-fraction variable (1997–2016) in the SEIB-DGVM SPITFIRE and GFED4s data, and most patterns differed only in the east region of Siberia (Figs. 4, S10). Compared to the burned-fraction variable, burned-area GFED4 has a smaller distribution pattern because it does not consider small fires (Fig. S9a). Comparison analysis of burned-fraction variables between SEIB-DGVM SPITFIRE and GFED4s showed a linear relationship with a correlation coefficient of R=0.87 (R²=0.75) (Fig. S11a). Similar to the comparison with GFED4s, the comparison of the SEIB-DGVM SPITFIRE output of burned-area variables with GFED4 data (1996–2016) shows a linear relationship with a correlation coefficient of R=0.78 (R²=0.61) (Fig. S11b). Furthermore, in the three regions (west, central, and east), the partial comparison of the burned-fraction variable with GFED4s showed values of R²=0.68, R²=0.51, and R²=0.58 (Fig. S13), while for the burned-area variable showed values of R²=0.51, R²=0.54, and R²=0.506 (Fig. S14), respectively. The burned fraction correlated better because both the GFED4s and the model's fire module considered small fires; scattered fire data with values less than 0.1 and approximately 0.1 were found in both the model's output and the GFED4s data.

The fire products (burned fraction and burned area) in the improved model have the same spatial distribution because they are calculated based on one core variable (fire probability) (Eq. 1). However, the spatial distributions of GFED4s (burned fraction) and GFED4 (burned area) differ for two reasons: first, because GFED4 does not consider small fires (Giglio et al., 2013), while GFED4s does, and, second, because GFED4s uses the modified burned-fraction equation, which is able to calculate the exact fire fraction and fuel load (not uniformized) in a grid cell (Van Der Werf et al., 2017).

Although the spatial distributions and patterns of the fire products (burned fraction and burned area) in the model and benchmark datasets (GFED4s and GFED4) data slightly differed, the model was able to produce annual mean value data that were similar to both benchmark datasets. The mean average burned fraction during 1997–2016 was 0.0137 in the simulations, compared to the GFED4s, which recorded the same value of 0.0137 with an RMSE value of $\mathrm{7.2}\times {\mathrm{10}}^{-\mathrm{4}}$. Furthermore, the mean average burned area of the model in 1996–2016 was 1428.5 ha per grid per year, compared to the GFED4 burned-area data, which closely recorded value of 1425.1 ha per grid per year by an RMSE value of 70.2 ha per grid per year. In summary, the model was able to produce mean average data that precisely resembled observational data.

Figure 4 Spatial distribution comparison of the annual average burned-fraction variable (1997–2016) from SEIB-DGVM SPITFIRE and GFED4s.

3.1.2 Aboveground biomass

The improved model simulated aboveground biomass values similar to those of the benchmark data. In 2010, 2017, and 2018, the simulations predicted 63.714±64.89, 64.141±65.54, and 64.313±65.61 Mg DM ha⁻¹ yr⁻¹, respectively, while the ESA Biomass CCI data showed 64.027±56.95, 64.548±54.69, and 65.05±55.78 Mg DM ha⁻¹ yr⁻¹, respectively, for the same years. The annual average AGB of the model in these years also showed the same increasing trend as the benchmark data, and the spatial distributions of the AGB model under CRU TS3.22 climate data and ESA Biomass CCI also agreed, with values of 83 %, 85 %, and 85 %, respectively (Figs. S15 and S16). Furthermore, when viewed on a smaller regional scale, the model is able to project better values in the west, central, and east regions, with average values of R²=0.73, R²=0.69, and R²=0.74, respectively (Fig. S17). Although there was an annual average increase in the number of forest fires, there was a high-variability trend in the model AGB values indicating succession after forest fires and a correct response to climate inputs variables based on each RCP scenario (Fig. 9d).

3.1.3 Annual and seasonal fluctuations in burned dry matter

The model's dry-matter data have a spatial distribution pattern similar to that of the model's fire products (burned fraction and burned biomass), as calculated from the available fire and fuel load data (fire product derivatives). The annual average dry-matter variability from the 1997–2016 model (under the historical climate product CRU TS3.22) and the GFED4s data agreed with 6.24 %, similar to the agreement of the fire products (Fig. S20). Spatial comparisons at the regional scale in the west, central, and east regions of Siberia show lower values than the Siberian region as a whole, with an agreement of 60.2 %, 64.4 %, and 58.8 % (Fig. S21).

We also compared seasonal dry-matter data to ensure that the monthly outputs of SEIB-DGVM SPITFIRE agree with the observations, as this difference influences seasonal aerosol emissions. Between 1997 and 2016, the GFED4s data exhibited high fluctuations and dynamics depending on the month and year, while SEIB-DGVM SPITFIRE was not able to reproduce these dynamics or accurately predict the occurrence of extreme events (Fig. S18a). For example, intense forest fires were recorded in 2003, 2012, and 2016. The monthly burned-dry-matter data for these years peaked in 2003 in May and in 2012 and 2016 in July (Fig. S18b–d). Severe wildfires in 2003 were due to low precipitation, as total precipitation reached only 36.0 mm in the Republic of Buryatia and 45.7 mm in the Chita Region between August 2002 and May 2003 (IFFN, 2003). The 41-year average precipitation between August and May (1981–2022) in the Republic of Buryatia was approximately 332.23 mm, and in the Chita Region was approximately 119.45 mm. Thus, the low precipitation in 2003 was an anomaly outside of the annual average range.

Furthermore, the improved model's monthly average burned dry matter in 2003, 2012, and 2016 was also lower compared to the GFED4s data. The burned-dry-matter values of the improved model were 58.64±5.86, 59.41±5.9, and 59.98±5.99 kg DM m⁻², while the benchmark data showed values of 122.36, 101.7, and 69.95 kg DM m⁻², respectively.

However, considering the entire period from 1997 to 2016, not only during years with extreme fire events but also for multiple years and time-series data was the model also able to reproduce similar average values. When comparing the monthly averages during 1997–2016, the model data yielded a value of 58.94±5.89 kg DM m⁻², while the GFED4s data yielded 59.12 kg DM m⁻². The model is not yet able to reproduce the exact value at a specific time of year or month because it runs in a long-term phase and is not yet able to predict sudden natural and anthropogenic conditions (factors). Overall, the spatial distribution comparison of the monthly dry-matter variables from GFED4s and SEIB-DGVM SPITFIRE for 20 years (1997–2016) revealed a correlation of 99 % (Fig. 5); therefore, the model was able to approximate the monthly averages.

3.1.4 Carbon dioxide (CO₂) and PM_2.5 emissions

Emissions from biomass burning contribute significantly to the global budget for residual gases and aerosols that affect the climate. It is estimated that biomass burning contributed up to 50 % of global CO and NO_x emissions in the troposphere (Galanter et al., 2000), and the most emitted gas during biomass burning is CO₂ (Ritchie et al., 2020). Since CO₂ emissions are the primary emissions that contribute to climate change, it is critical to assess and monitor them continuously.

In this study, out of 33 projected emissions (Tables 4 and S6), we validated the CO₂ variable that is able to represent all projected emissions because all estimated emissions are derived from the same burned-dry-matter variable, which differs only in the emission factor value of each gaseous emission. The highest annual average value of CO₂ emissions from 1997 to 2020 is from GFED4 data, followed by SEIB-DGVM SPITFIRE and then the GBEI product, with values of $\mathrm{105.64}\pm \mathrm{50.69}\times {\mathrm{10}}^{\mathrm{13}}$, $\mathrm{76.12}\pm \mathrm{0.87}\times {\mathrm{10}}^{\mathrm{13}}$, and $\mathrm{62.4}\pm \mathrm{26.09}\times {\mathrm{10}}^{\mathrm{13}}$ g CO₂, respectively (Table S3). The GFED4s and GBEI data have higher standard deviation values than the SEIB-DGVM SPITFIRE data and appear to have a large difference.

Spatially, the annual average CO₂ emission model data were 61.3 % (Fig. 6a) and 79.8 % (Fig. 6b) correlated with the GFED4s and GBEI data, respectively. Furthermore, CO₂ emissions of the model compared to the GFED4s in the three regions (west, central, and east) showed lower agreement than Siberia as a whole, at 62.7 %, 62.5 %, and 61.6 %, respectively (Fig. S29), whereas the comparison to GBEI data at the three regions showed agreements of 74.7 %, 77.6 %, and 64.3 %, respectively (Fig. S30). In addition, spatial comparison of annual mean data over 95 years (2006–2100) from SEIB-DGVM SPITFIRE, GFED4s, and GBEI datasets reveals similar values of 141.1±11.5, 157.2±14.8, and 148.7±7.12 Gg CO₂ yr⁻¹, respectively.

Our study area covers the boreal Asia (BOAS) area and a small part of central Asia (CEAS), differing from the GFED4 basis region classification; therefore, we extracted these areas from the GFED4s data for comparison (Fig. S22). A comparison of the GFED4s CO₂ data between the BOAS area and the Siberian area showed that the two datasets had a similarity of 98.2 % (Fig. S26), confirming the accuracy of the GFED4s validation data.

As all emission products are derived from fire products (dry-matter variables), emission factors displayed spatial and value dynamics similar to those of the fire products (Figs. S19, S27, and S43). When comparing the annual average dry-matter emission data and CO₂ emissions generated by the model, the results correlated perfectly (100 %, Fig. S31), indicating that the model runs well according to Eq. (3) and that the projected CO₂ and other emissions have the same distribution patterns as the dry-matter variable because all of the emission calculations are based on the dry-matter variable. However, they differ in their values because each emission species has a different emission factor.

We also compared the modeled PM_2.5 emissions and their distribution patterns with the Copernicus Atmosphere Monitoring Service (CAMS) (Romanov et al., 2022) data in seven Russian territories (Amur Region, Republic of Buryatia, Irkutsk Region, Khabarovsk Territory, Krasnoyarsk Territory, Transbaikal Territory, Yakutia (Republic of Sakha)) during 2010–2021. The improved model data and CAMS data both exhibited an increasing trend (Figs. 7a and 2 in Romanov et al., 2022) and a correlation of 85.8 % (Fig. 7b).

Figure 5 (a) Monthly average temporal-variation dry-matter emissions from GFED4s and SEIB-DGVM from 1997 to 2016. (b) Comparison of monthly average dry-matter emissions from GFED4s and SEIB-DGVM from 1997 to 2016. Standard deviation obtained from the monthly data from 1997 to 2016.

Figure 6 Latitude average spatial comparison of simulated CO₂ emissions from SEIB-DGVM SPITFIRE with (a) GFED4s from 1997 to 2016 and (b) GBEI from 2001 to 2020. Standard deviation obtained from the annual CO₂ emission data of each dataset.

Figure 7 (a) Temporal variation in projected PM_2.5 emissions under several climate scenarios from 2000 to 2100. Standard deviation obtained from the annual PM_2.5 emission value of each climate scenario (RCP8.5, RCP6.0, RCP4.5, and RCP2.6). (b) Comparison of PM_2.5 emissions from SEIB-DGVM SPITFIRE with the trend line of processed data from the Copernicus Atmosphere Monitoring Service in seven regions in Russia from 2010 to 2021 (Romanov et al., 2022).

3.2 Burned fraction

The improved model (SEIB-DGVM SPITFIRE) produces burned-fraction variables more accurately than the default model (SEIB-DGVM GlobFIRM). A spatial comparison of the average burned-fraction variables from 1997–2016 between GFED4s, SEIB-DGVM SPITFIRE, and the default SEIB-DGVM shows that SEIB-DGVM SPITFIRE achieves a 75 % similarity with GFED4s data, whereas the default model achieves only 68 % (Fig. S5).

The burned-fraction variable in the improved model exhibited a spatial distribution pattern different from that in the default model (Fig. S4a). According to the improved model, the burned-fraction data were distributed in the west, central, and southern areas (Fig. S4b). We compared the burned-fraction variable with the lightning flash rate and population density data to confirm that the produced variable considered the new ignition factor. The burned fraction showed a 46 % correlation with the lightning flash rate and a 6 % correlation with population density between 2006 and 2100 (Fig. S7a and b). In general, the burned fraction under all the RCP scenarios exhibited an increasing trend from 2006 to 2100, with the highest value occurring under the RCP4.5 scenario. Under the RCP4.5 scenario, the lowest value was 0.01371 and the highest value was 0.01427, with an average value of 0.01398 (Fig. S4d).

In contrast to the results produced from the improved model, the burned-fraction data from the default model were spread throughout most of the area (Fig. S4a). From 2006 to 2100 under all RCP scenarios, the burned fraction in the default model also exhibited an increasing trend. Under the RCP4.5 scenario, the lowest value is 0.002996 and the highest value is 0.003113, with an average value of 0.00306 (Fig. S4c), which is well below the outputs of the improved model.

3.3 Burned area

The burned area of the improved model showed a similar spatial distribution pattern under all the RCP scenarios (Fig. S6a). The distribution pattern of the burned-area variable was also similar to that of the burned-fraction variable, as the burned-area and burned-fraction calculation processes are both based on fire probability (Eq. 1). Overall, under all the scenarios, the burned area exhibited the same increasing trend, with the RCP4.5 scenario reaching the highest value. Under the RCP4.5 scenario from 2006 to 2100, the burned area has an average value of 1945.9 ha per grid per year and is projected to increase with values of 79.7 to 83.8×10⁵ ha (Fig. S6b). Since the default model does not compute burned area, this variable could not be compared between the improved model and the default model.

3.4 Burned biomass

The improved model confirmed uniform spatial distribution patterns for the fire variables: burned fraction (Fig. S4b), burned area (Fig. S6a), and burned biomass (Fig. 8b). All of the improved module fire variables confirmed to be mutually integrated because the calculation process comes from the first fire variable (burned fraction). Compared to the improved model, the spatial distribution pattern of the burned-biomass variable from the default model was wider and spread across the entire Siberian region (Fig. 8a). The spatial distribution pattern of the burned fraction (Fig. S4a) and burned biomass (Fig. 8a) in the default model is different and exhibited a box-like pattern in the center of the map. The internal model calculation flow relationship between the burned-fraction and burned-biomass variables in both the default and improved models shows a positive linear correlation, indicating harmony between these variables. A higher burned fraction corresponds to a higher burned biomass. The default model (SEIB-DGVM GlobFIRM) has an R² value of 0.83, while the improved model (SEIB-DGVM SPITFIRE) demonstrates better integration, with an R² value of 0.93 (Fig. S8a and d). Under all RCP scenarios from 2006 to 2100, the burned-biomass variable in both the default and improved models exhibited an increasing trend (Fig. 8c and d). This indicates correct integration between the burned-fraction and burned-area variables and an appropriate response to the climate input data. Furthermore, under the RCP6.0 climate scenario from 2000 to 2100, the burned-biomass value in the default model increases from 50.4 to 60.6 kg DM m⁻² (Fig. 8dc), while in the improved model it increases from 53 to 73.98 kg DM m⁻² (Fig. 8d). The 20-year variations and their trends of dry-matter emissions up to 2100 in the improved model (SEIB-DGVM SPITFIRE) are 55.90±1.31 (10.5 %), 60.52±1.12 (11.4 %), 64.43±1.36 (12.1 %), 69.23±1.37 (13 %), and 71.81±0.94 (13.5 %) (Fig. S32).

Figure 8 (a) Spatial distribution of annual average burned biomass from SEIB-DGVM GlobFIRM from 2000 to 2100. (b) Spatial distribution of annual average burned biomass from SEIB-DGVM SPITFIRE from 2000 to 2100. (c) Temporal variation in annual average burned biomass from SEIB-DGVM GlobFIRM from 2000 to 2100. (d) Temporal variation in annual average burned biomass from SEIB-DGVM SPITFIRE from 2000 to 2100.

3.5 Aboveground biomass

The aboveground biomass calculations in the default model and improved model used the same estimation process because the trunk biomass in SEIB-DGVM included coarse-root biomass; therefore, only approximately two-thirds of the trunk biomass was classified as aboveground biomass (Sato et al., 2007). However, during the calibration of the aboveground biomass variable with the ESA Biomass CCI benchmark dataset, we adjusted the calculation impact of fire and its distribution pattern (based on natural- and anthropogenic-ignition factors) on the availability of aboveground biomass.

According to the default model, the AGB distribution pattern appears to be the same as that of the fire variable; a box-like pattern still occurs on the map (Fig. 9a). Under the RCP8.5 scenario, from 2000 to 2100, the AGB increased from 63.72 to 120.1 Mg DM ha⁻¹ and the average value was 86.3 Mg DM ha⁻¹ (Fig. 9c). The aboveground biomass (AGB) variables in both the default and improved models exhibit an increasing trend and vary across RCP scenarios, with the highest values observed under RCP8.5 and the lowest under RCP2.6. This indicates that the models effectively read and process the RCP input climate data.

Compared to the default model, the improved AGB model has distribution patterns that are a bit different (Fig. 9b). In the central region of Siberia, some locations that have high AGB have been reduced due to the impact of forest fires so that the box-like pattern is no longer visible (Fig. 9b). The temporal variation in aboveground biomass in the improved model also shows an increasing trend due to the warming scenario of each RCP climate data input. The AGB under the RCP8.5 scenario from 2000 to 2100 increased from 59.08 to 126.7 Mg DM ha⁻¹ (Fig. 9d), and the mean was 88.68 Mg DM ha⁻¹. The 20-year variations and their trends of aboveground biomass up to 2100 are 65.45±1.19 (10.8 %), 71.69±2.90 (11.8 %), 83.38±3.61 (13.7 %), 99.17±5.06 (16.3 %), and 117.92±5.41 (19.4 %) (Fig. S33).

Figure 9 (a) Spatial distribution of annual average aboveground biomass from SEIB-DGVM GlobFIRM from 2000 to 2100. (b) Spatial distribution of annual average aboveground biomass from SEIB-DGVM SPITFIRE from 2000 to 2100. (c) Temporal variation in aboveground biomass from SEIB-DGVM GlobFIRM from 2000 to 2100. (d) Temporal variation in aboveground biomass from SEIB-DGVM SPITFIRE from 2000 to 2100.

3.6 Forest ecological variables under fire-on and fire-off simulations

We conducted complete simulations under fire-on and fire-off modes to compare and assess vegetation dynamics during forest fires. Assessing vegetation dynamics can be done by understanding the carbon pools in the certain region or globally, where carbon pools are easier to measure than carbon fluxes. In this study, the net primary production (NPP) is used as a reference variable because it is an important metric of the global carbon cycle (Running, 2022) and measures the rate of global plant growth. We obtained the NPP loss variable due to wildfire from fire-on and fire-off simulations. The NPP loss variable under all RCP scenarios shows an increasing trend. Under the RCP8.5 scenario, an average NPP loss of 385.19±40.4 g C m⁻² yr⁻¹ occurred during 2000–2100 (Fig. S25a). In addition to the NPP variable, the improved model (SEIB-DGVM SPITFIRE) can also simulate net biome production (NBP). Under the same RCP8.5 scenario, the annual average NBP from 2000–2100 shows a positive value of 307.7±43 Tg C yr⁻¹ (Fig. S25b), with a continuous increasing trend.

In relation to wildfires, assessing pre- and post-fire tree density variables is critical for measuring the impact of fires. Under the RCP8.5 scenario, in the fire-on simulation from 1997 to 2100, it is projected that the tree density in Siberia was 2181 trees ha⁻¹. However, under the same RCP and time range in the fire-off simulation, the tree density was 2363 trees ha⁻¹. We also compared the tree density between the fire-on and fire-off simulations under all the RCP scenarios and found that the tree density increased in the fire-off simulations compared to that in the fire-on simulations. Under the RCP8.5 scenario, on average, 174 trees ha⁻¹ yr⁻¹ died due to the fire (Fig. S25c).

We also conducted a more detailed assessment of several forest structure variables, such as the tree diameter at breast height (DBH), crown area, and tree height, from 2006 to 2100 under all the RCP scenarios. Under the RCP8.5 scenario, in the fire-on simulation, the results showed that tree DBH values varied from 0 to 4.7 m (average of 0.9 m), tree height varied from 0 to 75.4 m (average of 24.2 m), and crown area varied from 0 to 15.1 m² (average of 5.7 m²). The average tree structure in the fire-off simulation was greater than that in the fire-on simulation, with an average tree DBH, tree height, and crown area of 0.97 m, 24.1 m, and 6.5 m², respectively. The correlations between the tree structure variables under fire-on and fire-off simulation conditions were similar and highly correlated; the overall average correlation among the tree DBH, tree height, and crown area variables was 97 % (Fig. 10). Specifically, according to region classification, from the highest to the lowest, the values of tree height, tree DBH, and crown area are in the west region, central region, and east region. On average for 2081–2100 under RCP8.5 in each region, the tree height, tree DBH, and crown area variables show values of 28.43±0.8 m, 1.1±0.004 m, and 5.7±0.01 m² (west region); 28.3±0.9 m, 1.2±0.04 m, and 7.8±0.08 m² (central region); and 30.2±1.0 m, 1.2±0.06 m, and 8.5±0.2 m² (east region) (Figs. S37, S38, and S39). Furthermore, we found an interesting pattern; the simulated tree allometry variables (tree height, tree DBH, and crown area) in the east region of Siberia exhibit a greater range of values compared to those in the central and west regions of Siberia (Figs. S37, S38, and S39). Overall, all tree allometry variables in Siberia exhibit an increasing trend, and the differences between fire-on and fire-off simulations for all tree allometry variables are most pronounced in the east region of Siberia.

In addition, the relationship between the three variables (tree height, tree DBH, and crown area) in the west region and central region shows a linear trend where the higher the tree height, the greater the tree DBH and the wider the crown area (Fig. 11). The east region shows an interesting pattern, different from other regions, where there is low tree height (Fig. 11c). The west and central regions of Siberia exhibit a greater range of tree height values compared to the east region of Siberia (Fig. 11a and b). An interesting pattern was observed in the west region of Siberia, where trees with high tree height and a large DBH but low crown area were detected in some locations (Fig. 11a).

Figure 10 Annual average comparison of tree DBH, tree height, and crown area variables in (a–c) fire-on and (d–f) fire-off simulations (2006–2100). Each point represents one grid latitude average of each variable.

Figure 11 Relationships between simulated tree height, tree DBH, and crown area under the fire-on simulation and RCP8.5 scenario from 2000 to 2100 in the (a) west region, (b) central region, and (c) east region of Siberia.

3.7 Fire and AGB variable comparisons

We performed internal comparisons of fire and AGB variables within the improved model to ensure that the model worked properly and that the variable calculation processes were interrelated. The east region of Siberia had low fire patterns (Fig. S23), and, when compared with the AGB, this area also had very low AGB. We extracted the AGB data in the marked area with coordinates of 130–142° E, 65–80° N and discovered that the average simulated aboveground biomass in the area was 65.59 g C m⁻² from 1997 to 2023, compared to 416.4 g C m⁻² in the one-grid high-AGB areas. Furthermore, we assessed the fire danger index (FDI) variable in these low AGB areas and found that the mentioned region had a value of 0, indicating that it had a very low fire potential (Fig. S24a).

We also compared the fire variables (burned fraction, burned biomass) and AGB variables between the improved model and the default model. According to the default model, the correlation between the burned fraction and burned biomass was 0.83, the correlation between the burned fraction and AGB was 0.82, and the correlation between burned biomass and the AGB was 0.88 (Fig. S8a–c). According to the improved model, the correlation between the burned fraction and burned biomass was 0.93, the correlation between the burned fraction and AGB was 0.96, and the correlation between burned biomass and the AGB was 0.9 (Fig. S8d–f). Overall, both the default and improved models are well integrated, with the improved model demonstrating superior integration compared to the default model.

3.8 Future projection of burned-biomass emissions

Our model projects that from 2000 to 2100, regions in Siberia will produce CO₂ emissions ranging from 10 to 11 000×10⁸ g CO₂ yr⁻¹ (Fig. 12). The distribution patterns of CO₂ and other emissions are similar because all emissions are calculated based on the same variable dry-matter emissions. Over the 20-year period, we projected an increasing trend in CO₂ emissions across the various RCP scenarios, which aligns with the projected increase in forest fires through 2100.

The average from 2000 to 2100 shows that CO₂ emissions are highest under the RCP6.0, RCP2.6, RCP4.5, and RCP8.5 scenarios, with values of 885.8±75.4, 877.82±82.6, 871.4±80.6, and 865.5±69.6 Tg CO₂, respectively. Specifically, under the RCP6.0 scenario, the highest projected emissions are expected in the periods 2021–2040, 2041–2060, 2061–2080, and 2081–2100, with Siberia producing CO₂ emissions of 769.24±14.37, 830.52±15.61, 877.93±16.34, 940.46±20.59, and 981.73±12.61 Tg CO₂, respectively (Fig. 13).

The gaseous species with the highest emissions were CO₂, CO, PM_2.5, TPM, and TPC, and all of them exhibited similar increasing trends from 2000 to 2100 under all RCP scenarios. Under the RCP6.0 scenario, these emissions are expected to increase by 2.58±0.75, 0.21±0.06, 0.03±0.01, 0.02±0.01, and 0.014±0.006 Tg species yr⁻¹, respectively. The increasing trend of emission production until 2100 is also in line with the FDI variable, which shows the same increasing trend (Fig. S24b). Overall, by 2100, under the RCP6.0 scenario, the production of CO₂, CO, PM_2.5, TPM, and TPC emissions from forest biomass-burning combustion is projected to reach 1009.00±75.44, 80.74±6.04, 12.60±0.91, 9.54±0.69, and 6.61±0.48 Tg, respectively. The 20-year averages of the CO₂, CO, PM_2.5, TPM, and TPC emission data under all the RCP scenarios are provided in Table 4, and the other 28 emissions are provided in Table S6.

Figure 12 Spatial distribution of annual average projected CO₂ emissions (2000–2100) under the RCP8.5 scenario.

Figure 13 Temporal variation in projected CO₂ emissions from 2000 to 2100 under different RCPs scenarios.

Table 4 The 20-year average (±2 standard deviations) of projected emissions of CO₂, CO, PM_2.5, TPM, and TPC species from forest fires in Siberia (2023–2100). The emissions of the remaining 28 species are listed in Table S6. NA – not available.

4 Discussion

4.1 Feasibility of fire simulation

According to the default module, the fires spread throughout almost all of Siberia (Figs. S4a, 8a) because the module considered only the fuel amount and fuel moisture content. Thus, if the fuel load met the threshold requirement in any random grid, a fire appeared and could spread to other areas. Furthermore, the spatial distribution and trend of burned biomass under all of the RCP scenarios in the default fire module were not consistent with the burned-fraction data. Areas with high burned-fraction values should also have high values of burned biomass and vice versa.

However, in the improved module, the fires ignited only in areas that were covered in the lightning ignition and population ignition datasets based on the calculation of each ignition factor. This is confirmed by the comparison of the fire variable with the ignition factor variables. The comparison of the burned-fraction variable with the lightning flash strikes variable shows a strong correlation of 0.68 (R²=0.45), and the comparison of the burned-fraction variable with the population density variable shows a correlation of 0.24 (R²=0.06) (Fig. S7). These relatively low correlation values are due to the fact that the presence of an ignition factor does not guarantee that a fire will start; the area needs to have sufficient dry litter to feed the fire. Apart from these variables in the improved model, other factors also influence fire occurrence and spread in real life, such as slope and solar aspect (Rothermell, 1972), but their inclusion at this point was not possible due to the limitations of the model. In addition, when comparing the fire and AGB distributions, SEIB-DGVM SPITFIRE showed greater agreement than the default fire module.

However, differences remained between the spatial distribution patterns of the simulated fires and the GFED4s data in the east region of Siberia. We believe that the main reason for the lack of simulated fires in the east region of Siberia was the scarcity of available fuel and biomass for the ignition and spread of fires. We found that the AGB in these areas (130–142° E, 65–80° N; Fig. S23) was very low, averaging 65.59 g C m⁻². This value was far below the model minimum fuel load threshold requirement of 200 g C m⁻² (Sato et al., 2007) for fire ignition or spread. All three benchmark datasets, ESA Biomass CCI (aboveground biomass), GFED4 (burned area), and GFED4s (burned fraction), indicate that fire is present in this area, with ESA Biomass CCI showing an AGB of 2309.67 g DM m⁻². It is challenging to produce a model product that precisely predicts observations, as the simulations are highly dependent on the input data and dynamics, while the benchmark datasets were obtained from satellite image estimations that are able to capture natural conditions and events in real time. Even predictions based on satellite observations can differ significantly from field-based observations. For example, the International Forest Fire News (IFFN) for the Russian Federation reported that 2003 had extremely severe fires in Siberia based both on ground and aerial observation. However, the burned area was determined to be 2 654 000 ha based on field observations and 17 406 900 ha based on satellite-derived observations (NOAA AVHRR, Advanced Very High Resolution Radiometer) (IFFN, 2003; Siegert and Huang, 2005). The difference between ground observation data and satellite-derived data is due to differences in the data collection time and continuity. Ground-based observations are carried out only for a short time due to technical difficulties, while observations based on satellite data are carried out without any significant difficulties (IFFN, 2003). In this case, SEIB-DGVM SPITFIRE reported a burned area of 7 969 785 ha, an estimation centered between the observational and satellite data.

Overall, based on the fire variable outputs (burned fraction and burned area) from the improved model generated and validated with benchmark data, we project that Siberia will have an increasing trend until 2100 (Fig. S4d, S6b, and 8d). Yasunari et al. (2024) in a comprehensive assessment of the impacts of the Siberian wildfire using MIROC5 stated that there is high probability of increased Siberian wildfires in the future, and this estimate implies that worse air quality due to wildfires is predicted in the future, with frequent exceedances of air quality environmental standards (ESs).

Kasischke and Bruhwiler (2003) stated that the level of uncertainty in the burned-area parameter for estimating fire emissions in the Russian boreal forest is ±30 % for satellite imagery, while the uncertainty in the parameter is −300 % according to official government statistics, resulting in fires being largely underestimated. This difference in uncertainty was caused by the diverse parameters and equations used for estimation; the varying levels of detail of the analysis; and other factors, such as forest type, location, fuel load, fire type, and aboveground biomass density. Differences are also extrapolated when estimations for large areas are based on individual fires (Kasischke and Bruhwiler, 2003; Kukavskaya et al., 2013). Therefore, uncertainties will inevitably persist in model- or simulation-based research when comparing model- or simulation-based data with direct observations.

4.2 Forest resilience under fire and climate change

Terrestrial NPP is an essential element of the carbon cycle and global climate dynamics, as it directly affects the CO₂ content of the atmosphere, resulting in delayed climatic changes (Running, 2022). If NPP decreases, the land's ability to absorb CO₂ will decrease, causing atmospheric CO₂ to increase faster and thereby contributing to climate change (Running, 2022). Based on the comparison between the fire-on and fire-off simulation, under the RCP8.5 scenario, from 2000 to 2100 the NPP will decrease by 385.11±40.4 g C m⁻² (5.03±1.5 g C m⁻² yr⁻¹) due to wildfires until 2100. Satellite observations 1 year after boreal forest fires in Alaska and Canada recorded a 60–260 g C m⁻² loss of NPP (Hicke et al., 2003). In the coniferous forests of the western United States, post-fire NPP loss was also recorded and ranged from 67 to 312 g C m⁻² yr⁻¹ (Sparks et al., 2018). These data indicate that the NPP simulation results of SEIB-DGVM SPITFIRE are also consistent with some observational data in different areas.

NPP and NBP are both significant elements of the global C cycle and are used as indicators of ecosystem function and are linked to biodiversity, biogeochemical cycling, ecosystem resilience, and other aspects of ecosystem services (Pan et al., 2011; Richmond et al., 2007; Ito, 2011). However specifically, the mitigation ability of ecosystems is determined by net biome productivity (NBP) (Chapin et al., 2006; Fisher et al., 2014), and climate-driven large anomalies in NBP could impact the structure, composition, and function of terrestrial ecosystems (Frank et al., 2015). The 20-year average NBP from 2000–2100 shows a carbon sink in Siberia with an increasing trend (Fig. S34). Overall, from 2000 to 2100, RCP8.5 produces the highest value, followed by RCP6.0, RCP4.5, and RCP2.6, with values of 304.61±11.77, 286.78±10.64, 286.17±10.99, and 274.95±9.36 Tg C yr⁻¹, respectively (Table S4). The historical annual mean value of NBP in Siberia for 2000–2021 of 136.39±83.4 Tg C yr⁻¹ is also similar to the CLM4CN simulations (annual average of 1981–2006) in Eurasia and Boreal and Arctic of 204 and 284 Tg C yr⁻¹, respectively (Kantzas et al., 2013).

Under all climate scenarios from 2000 to 2100, we estimate that the net biome productivity (NBP) will continue to increase, indicating a continued flux of CO₂ from the atmosphere to the land (Fig. S34). The classification of NBP variables based on climate input data also shows the correct order, from the smallest under RCP2.6 to the largest under RCP8.5 (Fig. S35). This is because climate factors, such as temperature and precipitation, have a positive impact on vegetation (Yuan et al., 2021). On average, from 2000 to 2100, under the RCP8.5 climate scenario, the NBP in Siberia is estimated at 301.3±49.1 Tg C (equivalent to 3.01±0.5 Tg C yr⁻¹). Other studies have similar estimation that the NBP across northern peatlands, including the Russian Far East (RFE) and West Siberian Lowlands (WSL), ranges from 10 to 220 Tg C yr⁻¹ (Qiu et al., 2022). Additionally, we estimate that the heterotrophic respiration (HTR) in Siberia will continue to increase until 2100. On average under RCP8.5, from 2000 to 2100, the HTR value in Siberia is estimated at 4002.7±967.7 Tg C (equivalent to an increase of 40±9.7 Tg C yr⁻¹) (Fig. S36). We suggest that the high HTR values during those years were attributable to an elevated fuel load (Fig. S40a), followed by significant precipitation (Fig. S40b), which increased litter moisture content (Fig. S41a) and consequently accelerated the decomposition rates of litter and soil organic carbon. Increased heterotrophic respiration, tree mortality, and increased disturbance (drought and fire) contribute significantly to negative carbon fluxes from the ecosystem due to increased temperature and atmospheric CO₂ (Sharma et al., 2023). Overall, SEIB-DGVM SPITFIRE simulates that until the end of the 21st century, there will continue to be a strengthening of the land carbon sink in Siberia under all RCP scenarios. Boreal forests (1135 Mha) consistently acted as an average carbon sink of 0.5±0.1 Pg C yr⁻¹ over the 2 decades from 1990 to 2010. Furthermore, Asian Russia had the largest boreal carbon sink, which showed no overall change despite increased emissions from wildfire disturbances (Pan et al., 2011).

Boreal forest vegetation is naturally influenced by a variety of periodic disturbances, such as wildfires (Kasischke et al., 1995), insect outbreaks, and windthrow. Wildfires and insect outbreaks are not necessarily independent; there is a likelihood of wildfires often increasing or decreasing after insect outbreaks (Meigs et al., 2015, 2016). However, wildfires are among the main disturbances that drive forest dynamics, shape forest composition and structure, and affect biomass and productivity (Burns and Honkala, 1990; Greene and Johnson, 1999). Circumpolar northern boreal forests and tundra are likely to continue to warm more than most other terrestrial biomes according to available data from models and observations (Chapin et al., 2005; Foley, 2005; Meehl et al., 2007; Trenberth et al., 2007; Lee et al., 2021). Based on the observations and changes in regional attributes from 1950 to the present, it is projected that during 2071–2100, the WSB (west region of Siberia), ESB (east region of Siberia), and RFE (Russian Far East) will experience an increase in extreme temperatures with a high confidence of more than 7 °C for all seasons under the RCP8.5 scenario. Projected warming is most evident on the large continental Siberian Plateau, which has boreal and subboreal climates and biomes (i.e., taiga forests and tundra), during the winter season (Ozturk et al., 2017; IPCC, 2021). Such changes in climatic extreme scenarios and seasonality are also likely to have multiple effects, including extended but drier growing seasons, the occurrence of more intense convective storms leading to more lightning-caused fires (Hessilt et al., 2021; Kharuk et al., 2022), and decreased forest productivity (Orangeville et al., 2018); additionally, longer, warmer, and drier summers may cause an increase in fire frequency and size in some areas of boreal forests (Krawchuk et al., 2009; Flannigan et al., 2016; Wotton et al., 2017). This finding is in line with our results, which show that the assessment of forest ecology variables indicates tree mortality due to fire and succession as well as post-fire vegetation (Fig. S25c) affects NPP dynamics in Siberia (Fig. S25a).

Under the RCP2.6 scenario, SEIB-DGVM estimated the average tree density to be 2166 trees ha⁻¹ between 200 and 2023 in Siberia. The tree density is greater in northeastern Siberia (1197 trees ha⁻¹) than in southern Siberia (Miesner et al., 2022). Our simulation resulted in higher tree densities than did the observations in northeastern Siberia, as we covered a larger area of forest at 60–180° E, 45–80° N. The number of trees is affected by the frequency of fires at a certain location. Additionally, the number of trees destroyed by wildfires depended upon the climate scenario used in the simulations but naturally increased with fire frequency and size. In all the RCP scenarios, the number of destroyed trees was greater than that in the historical simulation and the number of destroyed trees increased annually, indicating that changes in climatic factors affected the density of surviving trees. The projected increase in the number of trees destroyed annually is consistent with the modeled fire product data, which exhibit an increasing trend until 2100. The difference in tree mortality data between climate scenarios is due to each climate scenario having a different projected temperature increase. In Siberia, under the RCP8.5 scenario, we simulate that the 2 m surface temperature will increase by 4.67 °C by 2100 (Fig. S42). This estimate aligns with the IPCC projections, which predict air temperature increases by 2100 ranging from 0.3–1.7 °C (average of 1.0 °C) under the RCP2.6 scenario, 1.1–2.6 °C (average of 1.8 °C) under the RCP4.5 scenario, 1.4–3.1 °C (average of 2.2 °C) under the RCP6.0 scenario, and 2.6–4.8 °C (average of 3.7 °C) under the RCP8.5 scenario (IPCC, 2021).

The DBH ranges of the trees in the fire-on and fire-off simulations were comparable to those in northeastern Siberia, where the DBH ranged up to 71.6 cm, the tree height ranged up to 28.5 m, and the crown area averaged 4.77 m² (Miesner et al., 2022). As the average DBH variable was similar in the fire-on and fire-off simulations, trees with large DBHs are resistant to fire. This was also confirmed based on observational research in Yenisei Siberia, where trees with a DBH greater than 18.1 cm were the most resistant to further post-fire succession (Bryukhanov et al., 2018). Specifically, based on the division of regions, we found an interesting pattern, in which the east region of Siberia has the highest value of allometry variables (tree height, tree DBH, and crown area), followed by the central region, with the lowest in the west region (Figs. S37, S38, and S39). An interesting pattern was observed in the west region of Siberia, where trees with high height and a large DBH but low crown area were detected in some locations (Fig. 11a). We suggest that this happens because of the wildfire; the central region of Siberia has the highest wildfire frequency, followed by the west region, then the east region. The major Siberian forest types are formed by larch (Larix sibirica, L. gmelinii, and L. cajanderi), with the majority distributed in the west and central regions of Siberia (Fig. 1, and Fig. 1: Kharuk et al., 2021). Furthermore, larch is classified as pyrophytic species that have adapted or evolved under conditions of periodic forest fires; they have adapted and gained a competitive edge over non-fire-adapted species in regenerating and growing in burned areas (Kharuk et al., 2021). The abundance of species and high frequency of wildfires in the central and west regions of Siberia led to excellent larch succession and regeneration as evidenced by the high tree allometry variables, but the projected continuity of wildfires led to a downward trend. On the other hand, in the east region, very few wildfires are simulated, partly due to the low aboveground biomass available in some areas, which affects ignition and fire spread (Fig. S23). However, due to the low frequency of wildfires, allometric variables are projected to have an increasing trend until 2100 in the east region of Siberia. The unique relationship between allometric variables, which are naturally distributed without a wide gap between grid plots (Fig. 11c), in the east region is also due to the area's low wildfire frequency. The forest in the east region of Siberia appears to grow and spread naturally without a high impact from the wildfire. While the majority of tree species in Siberia (larch) regenerate extremely well on post-fire-mineralized soil, they regenerate very slowly over a ground floor covered in lichen and moss (where the soil's surface is tough for sapling roots to reach) (Kharuk et al., 2016).

4.3 Spatial distribution and temporal variation in biomass-burning emissions under climate change scenarios

The spatiotemporal dynamics of the biomass-burning emissions under all RCP scenarios had similar patterns and trends, but they had slightly different variations in dynamics because climate affects the frequency and distribution of fires. This is evidenced by all fire variables produced by the model, from the burned fraction to burned-biomass emissions. In the last 20 years of the projection (2080–2100), the highest values were obtained from simulations using the climate inputs of RCP2.6, RCP4.5, RCP6.0, and RCP8.5. This occurs because each RCP scenario exhibits varying radiative forcing, with RCP8.5 notably experiencing the highest temperature increase (Fig. S42) and also projecting the highest precipitation levels (Fig. S40b). The fuel load variable follows a corresponding order reflective of RCP forcing levels, with RCP8.5 showing the highest and RCP2.6 showing the lowest (Fig. S40a). However, due to increased precipitation and temperature-induced snowmelt, the moisture content of litter fractions in RCP8.5 simulations attains the highest values, contrasting with the lowest values in RCP2.6. Consequently, available fuel loads may not ignite in areas with high moisture content, leading to projections of the highest burned-biomass emissions in the last 20 years of RCP climate projections (2080–2100) for RCP2.6, RCP4.5, RCP6.0, and RCP8.5, respectively. The difference in emission values between climate scenarios in the same year shows that temperature has an impact on vegetation succession and climate-sensitive emission production from wildfires (Gutierrez et al., 2021; Stocker et al., 2021). Thus, the model is able to simulate and integrate fire disturbance, forest dynamics or vegetation succession, and burned-biomass emissions well.

Over a 20-year average from 2080 to 2100, under RCP6.0 in Siberia, our simulation predicts that forest fires will emit CO₂, CO, PM_2.5, TPM, and TPC in amounts of 989.93±13.72, 79.24±1.10, 12.91±0.18, 9.78±0.14, and 6.77±0.09 Tg, respectively (Table 4). Spatially, the projections depict heterogeneous patterns of burned-biomass emissions, with regions of high emission intensity concentrated in areas of larch forest (Larix spp.), consistent with Fig. 1 and our simulation results, where the fire and emission variables show high values in the central region to southern Siberia (Figs. S4b, 8b, and S6b). This is reinforced by field-based estimation data showing that fires in this region result in high tree mortality (76 %), with Siberian larch forests experiencing greater aboveground carbon loss after fire than North American forests, both in absolute and relative levels (Webb et al., 2024). We also visualized all the 33 graphs depicting projected burned-biomass emissions, offering valuable insights into the future dynamics of the burned-biomass emissions in Siberia. Across these graphs, we observe distinct temporal patterns, revealing trends in burned-biomass emissions over time. Under the RCP8.5 to RCP2.6 scenarios, the 20-year average comparison of overall burned-biomass emission data from 2080–2100, compared to data from 2000–2020, shows projected increases of 23.87 %, 27.63 %, 29.34 %, and 30.36 %, respectively (Fig. S43). The 20-year dynamics are summarized in Tables 4 and S6. Furthermore, each year, various climate scenarios predicted differing emissions based on the respective radiative forcing values (from lowest to highest). The RCP8.5, RCP6.0, RCP4.0, and RCP2.6 scenarios exhibited average annual increases of 0.295 %, 0.354 %, 0.358 %, and 0.361 % yr⁻¹, respectively, from 2000 to 2100.

Under the RCP4.5 scenario, radiative forcing stabilized until 2100 (Thomson et al., 2011), which is consistent with our results, as emissions under the RCP4.5 scenario were more stable than those under the other RCP scenarios. Therefore, it is indicated that the trend in fire emissions is consistent with the different scenario-dependent trends in radiative forcings. Overall, based on the RCPs climate scenario data used (MirocAR5), the emission scenario projected an increase in global mean surface temperature in the range of 1.0–3.7 (0.3–4.8) °C (IPCC, 2014) and currently ranges between 1.5 and 6.0 °C by 2100 compared to 1850–1900 mean value (Lee et al., 2021). One of the impacts of rising global temperatures is the increased occurrence and severity of forest fires, which lead to a greater prevalence of wildfire (Schoennagel et al., 2017; Haider et al., 2019). The global land area burned by wildfires is expected to increase by 35 % if the global temperature increases by 2 °C and precipitation patterns change (Pörtner et al., 2022). Extremely high temperatures increase the frequency of severe droughts and proliferate wildfires in several regions, such as southern Europe, northern Eurasia, the USA, and Australia (IPCC, 2021). These frequent and severe wildfires will inevitably lead to an increase in the atmospheric concentration of biomass-burning products (Marlon et al., 2008; Amiro et al., 2009; Tian et al., 2023).

Forests in Siberia are very important for continuous monitoring and assessing because they have a significant impact on regional (short-term) and global (long-term) air quality and human health due to the large amounts of carbon emissions, smoke aerosols, and trace gases in the atmosphere. In addition to the observed amount of emissions, organic carbon–elemental carbon emissions exceeded 3 times greater and emissions of inorganic ions (SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and NH${}_{\mathrm{4}}^{+}$) were found to be 5 times greater than the annual average wildfire emissions from August 2010 to August 2011 (Popovicheva et al., 2014). Increased Siberian wildfire aerosols would significantly degrade air quality, particularly in the surrounding and downwind regions of Siberia (Yasunari et al., 2024). The emitted substances can be transported over long distances and affect air quality in other regions, including North America and northeastern China (Teakles et al., 2017; Johnson et al., 2021; Sun et al., 2023).

Estimating future fire emissions and their impact on air quality is challenging due to model limitations and uncertainties in estimation methods, potential mixing of emissions in the atmosphere, climate radiative forcing factors, and emission transport (Winiger et al., 2017; Schacht et al., 2019). SEIB-DGVM SPITFIRE was not able to reproduce the events in the validation data for the same year or month but simulated similar dynamic patterns and values. This difference occurs because the benchmark data obtained from satellite image data closely follow natural conditions, while the model accumulates uncertainties due to its long simulation period. The emission estimation method used in the model refers to the dry-matter variable and the emission factor from Andreae and Merlet (2001) and Andreae (2019), where the emission factors are obtained from laboratory and small field experiments. Each species has specific characteristics that require different assessment methods, and the combustion characteristics can be very different from those of large-scale open biomass burning and wildfires. Kasischke and Bruhwiler (2003) reported that the level of uncertainty in the emission factor parameters for estimating emissions from fires in Russian boreal forests was ±20 %–50 %, which agrees with the ±50 % uncertainty level for major emissions presented by Andreae and Merlet (2001). However, in SEIB-DGVM SPITFIRE, we also used the latest emission factor from Andreae (2019), which was developed for oxygenated volatile organic compounds and for HCN; this approach improved all assessment compound emissions significantly with more accurate measurements and has been widely used by various dynamic global vegetation models to estimate biomass-burning emissions globally. Overall, the comparison between different climate RCP scenarios provides further insight into uncertainties and variability in the projections, offering valuable information for understanding the potential impacts of future burned-biomass emissions on air quality, climate dynamics, and ecosystem health. Through this analysis, our study contributes to a better understanding of the drivers and implications of burned-biomass emissions, informing policy decisions and management strategies aimed at mitigating their environmental and societal impacts.

4.4 Model uncertainty

Our study is a process of combining and improving SEIB-DGVM with the SPITFIRE fire module, each of which has different characteristics and some default variables, parameters, and inputs. The implementation of any complex model improvement inherently introduces uncertainty, stemming from various sources such as parameterization choices, model structure, new model input, and the representation of complex biophysical processes. Specifically in the context of the fire module enhancement, uncertainties may arise from the characterization of ignition sources, fire behavior, fuel dynamics, and fire spread mechanisms. These uncertainties can significantly influence the accuracy and reliability of model projections, particularly in simulating the spatial and temporal patterns of fire occurrence, intensity, and impacts on vegetation dynamics and carbon cycling.

Kasischke and Bruhwiler (2003) mentioned that, in some cases, the data needed to generate input parameters for those equations are very well defined, whereas in others, they are based on a very limited set of the observations data. Thus, the input data selection and the parameter setting for those equation calculations are the source of the uncertainty, with the provided emission range (Table 5 in Kasischke and Bruhwiler, 2003). In more detail Kasischke and Bruhwiler (2003) stated that uncertainty can be classified into two groups. First, there are environmental characteristics (direct or indirect observations) including location, fire type, and aboveground/ground layer biomass. Second, there are uncertainties from parameters that can be measured in individual biomass combustion processes, while the application on a large scale, with time differences and climatic influences, is very challenging, and the combustion process consists of several stages, with different emissions for each stage.

To mitigate model uncertainty, we have employed rigorous model verification, calibration, and evaluation procedures, comparing model outputs with several benchmark datasets (Table 3). The verification helps to ensure the new inputs (ignition factors) can be read, processed, and output properly. The process of calibrating all major variables with benchmark datasets is carried out sequentially and with several iterations, ensuring that the output of individual variables matches the benchmark dataset that is the target of validation. This validation process helps assess the model's ability to reproduce historical fire patterns and dynamics accurately.

However, this does not eliminate the uncertainty in the model we have developed, and we still have limitations where the emission variable distribution pattern of emissions is strongly influenced by the pattern of the resulting fire because the emission variable is calculated with the same dry-matter emissions. This spatial distribution also affected other vegetation variables, due to the relationship calculation between fire and vegetation variables. This is a potential further study for adjusting the factors that affect the distribution of fire to be similar to the benchmark data. In addition, natural factors that affect the dynamics of significant fire disturbances at specific times are still not well simulated by our model. Inversely, our model also has the advantage of being able to numerically simulate averaged data in the long term very accurately; based on the results of numerical comparisons with benchmark data, the model is able to simulate with a value of 99 %.

5 Conclusions and recommendations

We introduced the SPITFIRE fire module into SEIB-DGVM and achieved a better representation of fire dynamics in Siberia between 1996 and 2100 by creating monthly outputs and producing several new outputs related to fires at a 0.5° spatial resolution, such as vegetation and burned-biomass emission variables. Our modifications have led to a more realistic depiction of fire frequency, intensity, and extent, aligning the model outputs more closely with benchmark datasets. Compared to the default model, the improved model (SEIB-DGVM SPITFIRE) demonstrates a higher accuracy in simulating the burned fraction, achieving a 75 % agreement with the GFED4s data, whereas the default model only shows a 68 % agreement. Overall, the major fire-related variables (vegetation, CO₂ and PM_2.5 emissions, burned area, burned fraction, aboveground biomass, and dry matter) all achieved an average spatial agreement of 70.7 % or higher with the observational data. Additionally, the improved model accurately simulated forest structure, increasing the agreement between the simulated and observed dataset patterns and further emphasizing the reliability of our model and its emission projections. Under the RCP2.6 scenario, we estimated that the CO₂, CO, PM_2.5, total particulate matter (TPM), and total particulate carbon (TPC) emissions in Siberia will continue to increase annually until 2100 by an average of 2.71±0.87, 0.22±0.07, 0.04±0.01, 0.03±0.01, and 0.02±0.01 Tg species yr⁻¹, respectively. Moreover, forest fires in Siberia in 2100 are projected to emit all five of these compounds under the RCP8.5 scenario, amounting to 1010.00±82.64, 80.84±6.61, 13.17±1.08, 9.97±0.82, and 6.91±0.57 Tg, respectively.

Although our research has taken significant steps, there are several limitations that require further research. Future studies should minimize the uncertainty in the simulations and achieve better fits with benchmark datasets of fire, vegetation, and emission products. Specific parameter settings need to also be developed to emphasize regional and seasonal differences. Continued improvement in the fire module and consideration of feedback loops will be crucial for continuously enhancing the accuracy of our models. Our work contributes to a more comprehensive understanding of the intricate interactions between fire dynamics, ecosystems, and climate, creating a new path for informed decision-making and broadening the field of biogeochemistry and global elemental cycles, as well as highlighting the importance of accurate vegetation dynamic modeling.

Appendix A: Input and output of SEIB-DGVM SPITFIRE

A1 Input

The following is the input of SEIB-DGVM SPITFIRE:

1. location – latitude and altitude;

2. soil (fixed in time) – soil moisture at the saturation point, field capacity, matrix potential, wilting point, and albedo;

3. climatic data (daily) – air temperature, soil temperature, fraction of cloud cover, precipitation, humidity, and wind velocity;

4. atmospheric carbon dioxide (CO₂) concentrations;

5. fire ignition factors – population density (GPWv4) and lightning flash rate (LIS/OTD HRFC).

A2 Output

The following is the output of SEIB-DGVM SPITFIRE:

1. carbon dynamics (daily–yearly) – terrestrial carbon pool (woody biomass, grass biomass, litter, and soil organic matter) and CO₂ absorption and emission rates;

2. water dynamics (daily) – soil moisture content (in three layers), interception rate, evaporation rate, transpiration rate, interception rate, and runoff rate;

3. radiation (daily) – albedo from the terrestrial surface;

4. properties of vegetation (daily–yearly) – vegetation type, dominant plant functional type, leaf area index, tree density, size distribution of trees, age distribution of trees, woody biomass for each tree, and grass biomass per unit area;

5. disturbances (monthly–yearly) – fire fraction, burned area, burned biomass, FDI, complete SPITFIRE variables, and 33 types of burned-biomass emissions.

Table A1 Processes in SEIB-DGVM SPITFIRE and the approaches used to represent each process.

Table A2 PFTs, variables, parameters, and constants in the model's equations.

Additional equations and variables of the implemented SPITFIRE module are adjustments to Thonicke et al. (2010) and are found in Table A2 and the Supplement.