A bottom–up savanna fire fuel consumption inventory and its application to savanna burning in Kafue National Park, Zambia

Case study area

Kafue National Park (15°6′4.0″S, 25°41′20.7″E) lies in western Zambia, covering an area of approximately 22,400 km² across Southern, Central and North-Western provinces. It is one of the largest parks or protected areas in Africa (UNEP-WCMC and IUCN 2023), comparable in size with smaller countries such as Rwanda or Slovenia. The park covers a rainfall gradient ranging from as low as 500 mm year⁻¹ in extreme cases at the southern end to upwards of 1100 mm year⁻¹ further north (Hersbach et al. 2020). The region has two distinct seasons, the wet season (November–April) in which the bulk of the annual rainfall occurs and the dry season (May–October) when most biomass burning occurs. The park is surrounded to the north, east and south by Game Management Areas (GMAs) and some forest reserves. The Lusaka–Mongu highway bisects the park longitudinally, and the Kafue river flows along much of the eastern edge and up through the northern region of the park. Land cover consists primarily of a combination of miombo woodland and (occasionally flooded) grasslands (dambos), with some riverine areas close to (seasonal) riverbanks (Olson et al. 2001). Dominant miombo tree species include the legumes Julbernardia paniculata, Brachystegia spiciformis and Burkea africana. The park is home to over half of Zambia’s bird species, dozens of fish and amphibian species live in its waterways, and over 150 different species of mammal including lions, elephants, cheetah and wild dogs (Mkanda et al. 2018) can be found within its borders.

Emissions calculations overview

Biomass burning emissions are calculated using the following equation from Seiler and Crutzen (1980):

where emissions are in grams per species, FL is fuel load in kilograms dry matter per square metre, BA is burned area in square metres, CC is the fractional combustion completeness of the fuel, and EF is the emission factor in grams per kilogram dry matter (DM). We built the FL inventory from the bottom up using a library of in situ field measurements. Each field site, containing two or three transects, was then surveyed with UAV-mounted cameras. The field sites were distributed across a rainfall gradient (Fig. 1) but also capturing some different soil types. Care was taken to place sites away from land used for active grazing (e.g. cattle posts or other animal husbandry), and the effects of ambient wild grazer activity are by nature included in measurements. Field sampling procedure details are given in Russell-Smith et al. (2021). A number of data layers were extracted from each site survey, and fed into a random forest (RF) regressor to create UAV-generated FL maps on a transect scale. These maps were then gridded onto a 20 m grid, each grid cell matched with a Sentinel-2 pixel, and these data were run through a second RF regressor to create Sentinel-2 scale FL maps. This process is described in the Fuel load sub-section. We compared emissions generated using several different BA datasets in this study – a Sentinel-2 based product (Chuvieco et al. 2022), the long-term MODerate Resolution Imaging Spectrometer (MODIS)-derived product (Giglio et al. 2018) and Landsat-8 based manually constructed dataset specifically for KNP. This last dataset was built by thresholding difference indices, and then filtering the resulting layers for anomalies. These datasets are addressed in more detail in the Burned area sub-section. The methods and data used to create FL, BA and CC layers are detailed in the following sections. EF data were generated on a global grid scale using extensive smoke sampling from a UAV system, augmented with other measurements from throughout the literature. Coupled with vegetation and meteorological reanalysis data, monthly dynamic EFs on a 0.25° × 0.25° grid were produced (Vernooij et al. 2023).

Fig. 1.

Distribution of plots in southern Africa across study sites in Tsodilo Hills (NW Botswana), Kasane Forest Extension (NE Botswana), Bovu Protected Forest Area (SW Zambia), Kafue National Park (dotted line, W Zambia), Lualaba Protected Forest Area (NW Zambia) and Niassa Special Reserve (N Mozambique). Early dry season plots have blue markers, and late dry season plots have red markers. Diamond-shaped markers represent plots established in 2019, triangular-shaped markers represent plots established in 2021 and circular markers those established in 2022.

Fuel load

The process of creating detailed biomass maps for KNP began with building an inventory through in situ measurements. This consisted of measuring different categories of fuel within groups of transects in a single plot, which was simultaneously surveyed with a UAV-mounted camera. Orthophotos were created encompassing all transects within a plot, and individual transects were demarcated within this orthophoto. Using UAV-based features and the meteorological history of the plot, an RF model was created to predict biomass for the entire orthophoto (hereafter UAV model). These biomass orthophotos were then distributed across Sentinel-2 pixels to form a training set for a second RF model, which predicted biomass on the basis of exclusively Sentinel-2-derived features (hereafter S2 model). This was then used to create a map of monthly biomass in KNP, on the basis of Sentinel-2 imagery. These steps are detailed in the following sections.

UAV model

Two UAV systems were used to survey the plots. The first was a DJI Matrice 100 with a mounted MicaSense RedEdge Multispectral camera, used to survey all plots laid in 2019. The second was a DJI Phantom P4 Multispectral, used to survey 2021 and 2022 plots. Both systems made use of the same software for image capture and processing (DJI Groundstation Pro and DJI Terra respectively) and the flight parameters were kept the same for both systems: altitude was fixed at 50 m and image overlap set to 80%, such that after stitching, the orthophotos had a resolution of 5 cm. Flight times were restricted to between 10:00 and 15:00 hours daily so as to avoid flying under low sun angles. Weather conditions while sampling were almost exclusively cloudless. Both cameras share the same five bands – blue, green, red, red edge and near infrared. In some bands, there is less overlap between the two cameras than in others, although for the longer wavelengths predominantly used in this study, the effect is less pronounced. Although we do expect there to be some difference in vegetation indices (VIs) calculated from bands that are not identical (Teillet et al. 1997), there is generally a strong correlation between VIs across different UAV-mounted instruments (Deng et al. 2018) and the measurements were taken in regions with similar features and under similar weather conditions. Furthermore, much of the UAV-generated data we rely on in this study are band-independent. Classification maps (Fig. 2c) and digital elevation models (DEMs, Fig. 2d) were generated individually for each orthophoto using the same algorithms as a method of fully unifying data captured by two similar but separate systems. For FL quantification, fuel class maps were generated using Geographical Object-Based Image Analysis (GEOBIA, Clewley et al. 2014) and combined with selected VIs and meteorological history from ERA5 Land hourly data (Hersbach et al. 2020). Each transect was ascribed a set of features listed in Table 3. Descriptions of the processes involved in obtaining these features are found in Eames et al. (2021). Examples of features obtained from UAV systems after post-processing are shown in Fig. 2. Once all features for each transect had been collected, they were fed together with the fuel measurements for each fuel category (Table 1) into an RF model (Breiman 2001; Geurts et al. 2006) implemented in Python (Pedregosa et al. 2011) with 10-fold cross-validation. Hyper-parameter tuning was done using a randomised grid method. A total of 20% of data points was kept aside for validation purposes. Tree biomass was developed separately by directly linking SAR backscattering from the Advanced Land Observing Satellite Phased Array type L-band Synthetic Aperture Radar (ALOS-PALSAR) to tree biomass, in combination with LiDAR validation and ground-based measurements (Bouvet et al. 2018). The degree of RADAR backscattering is directly linked to the amount of standing vegetation in dry and open biomes such as savannas. Bouvet et al. (2018) matched in situ measurements from plots across the African continent with backscatter data, and related them via a semi-empirical exponential equation corresponding to the modified Water Cloud Model (Santoro et al. 2002). The 25 m source data were regridded to 20 m using a weighted averaging algorithm within the boundaries of KNP to match with the Sentinel-2 grid. Following on from the UAV model creation, a 20 m Sentinel-2 grid was overlaid on each orthophoto, and features in Table 3 were extracted and calculated for each grid cell fully within the orthophoto. The UAV model was then applied to these grid cells, resulting in a 20 m resolution biomass map. An example of this is shown in Fig. 2e. Each grid cell was paired with a Sentinel-2 pixel from the most recent overpass previous to the UAV flight. The Sentinel-2 twin satellites provide a 5-day revisit time over Southern Africa, though in some rare cases, cloud cover and cloud shadow resulted in the most recent image being unusable. In these cases, the next available cloud-free image was used, up to a maximum period of 1 month before the UAV flight date. The total fine fuel class performed approximately as well in the UAV RF model as grass and better than litter (shown in Results section). Consequently, only the total fine fuel class was used to train the S2 model instead of using separate models for grass and litter. However, litter will continue to be included in its function as a proxy for coarse woody debris.

Fig. 2.

Example orthophoto from the Phantom P4 Multispectral camera in Kafue National Park (15°24′26.6″S, 25°55′12.5″E). The red-green-blue (RGB) image is shown in (a) along with an example transect (red box). The Normalised Difference Vegetation Index (NDVI) for this orthophoto is shown in (b), the GEOBIA classification in (c), and DEM in (d). Red areas in the DEM are equivalent to an elevation of 0 m (ground). Panel (e) shows the UAV model output for grass on a Sentinel-2 20 m grid.

Table 3.Features used in building the UAV RF model. ‘Last fire’ refers to the date of the most recent fire detected in MODIS MCD64A1 C6.1 data.

Feature name	Source	Description
Digital Elevation Model (DEM)	UAV	Mean canopy height from a photogrammetry-derived elevation model
Burned Area Index (BAI)	UAV	Mean burned area index value (Martín et al. 2006)
Leaf Chlorophyll Index (LCI)	UAV	Mean leaf chlorophyll index value (Datt 1999)
Normalised Difference Vegetation Index (NDVI)	UAV	Mean NDVI value (Tucker 1979)
Near infra-red (NIR)	UAV	Mean reflectance in the NIR band
Proportion of bare soil	UAV	Fraction of the transect classified as bare soil
Proportion of foliage	UAV	Fraction of the transect classified as foliage of live leaf material
Proportion of grass	UAV	Fraction of the transect classified as grass-covered
Proportion of wood	UAV	Fraction of the transect classified as woody material, live or dead
Sum of precipitation	ERA5	Total precipitation (mm) since the last fire
Mean annual precipitation	ERA5	Mean annual precipitation (mm) in the period since the last fire
Mean surface solar radiation	ERA5	Mean solar radiation (J m−2) incident on the transect since the last fire
Total evaporation	ERA5	Total depth (m) of water equivalent evaporated from the transect since the last fire
Shrub fraction	UAV	Fraction of the transect classified as shrub
Shrub percentage	UAV	Fraction of the total canopy within the transect classified as shrub

Burned area

Three different sources of BA data were used for the study area in 2019. The lowest-resolution dataset used was the 500 m MODIS MCD64A1 Collection 6.1 dataset (Giglio et al. 2018), available on a daily basis. The second was FireCCISFD20 (Chuvieco et al. 2022), a 20 m product based on images from the twin Sentinel-2 satellites. The final BA dataset (hereafter referred to as the Yates dataset) was produced manually for 90% of KNP’s area from 30 m Landsat imagery using methods similar to the North Australia and Rangelands Fire Information project (NAFI, https://firenorth.org.au/). Using indices incorporating the NIR and SWIR sections of the spectrum that are sensitive to changes in vegetation such that burn scars are easily detectable (Garcia and Caselles 1991; Pereira et al. 1999), these potential burn scars were mapped using a threshold, and then examined visually to identify and filter areas where the indices were confounded by non-fire related processes, e.g. surface water changes including flooding or drought (Evans and Russell-Smith 2020). The narrow western section of the park not covered by the Yates dataset was supplemented by FireCCISFD20 BA. As this dataset was developed under scrutiny of the KNP area in particular, it will be assumed to be the most accurate for the remainder of this paper. A comparison of the different BA products listed, including a MODIS-based BA with statistically modelled ‘small fires’ from GFED5 (Chen et al. 2023), is shown in Fig. 3.

Fig. 3.

Burned area from satellite-derived products in the year 2019. The red line is FireCCISFD20 derived from the Sentinel-2 MSI, the blue line is MODIS MCD64A1 C6.1, the yellow line MODIS-based GFED5 burned area including small fires, and the black line the manually compiled Yates dataset based on Landsat-8 imagery. The blue dotted line and corresponding right-hand side y-axis shows the mean monthly rainfall for KNP in 2019 (Hersbach et al. 2020). In the peak burn month (June), approximately one-quarter of the park area is burned.

Combustion completeness

Combustion completeness (CC) is the fraction of available FL consumed in the process of biomass burning. This varies according to fuel class, fuel/soil moisture, and weather conditions (van Leeuwen et al. 2014; Russell-Smith et al. 2021). Generally speaking, surface fuels burn most completely, and larger surface-level objects such as coarse woody debris or young shrubs are consumed to a lesser degree, either because of retained moisture in live standing fuel or because fires pass through too quickly to ignite larger objects. Mature trees are consumed the least and usually remain unaffected by less severe fires (Russell-Smith et al. 2009). It is difficult to account for local weather and wind conditions that might affect CC for one particular fire, as these can be very small scale and change rapidly during the course of the fire itself. This makes mapping CC on a pixel level exceedingly challenging. CC is, however, related to environmental water content (Shea et al. 1996 – see also Supplementary Material S4, Fig. S6). We therefore use SM data again as a basis for a general CC map for KNP. We broadly follow the methods used in van Wees et al. (2022), employing an SM-derived scalar (ϵ_SM, Eqn 2) to scale CC. In our case, the SMAP data were averaged to monthly values to scale between a CC range for each fuel class (Table 4). These ranges were adjusted from van Wees et al. (2022) and augmented with other field data from van Leeuwen et al. (2014) to match CC ranges from field measurements within KNP. We acknowledge that this does not account for all the factors influencing CC, but in this context represents an adequate solution to the challenge continuous CC representation in large-scale modelling presents. The scalar used for pixel x and time t is:

where 0.45 m³ m⁻³ is the typical maximum SM outside wetland areas (van Wees et al. 2022). The role of fire in relation to trees in tropical savannas is almost always as a sapling removal mechanism rather than directly affecting adult trees (Higgins et al. 2000; Ojeda et al. 2005) and in KNP in 2019, Global Forest Change data (https://glad.earthengine.app/view/global-forest-change) based on 30 m Landsat imagery showed little forest loss, and none due to fire (Tyukavina et al. 2022). We therefore expect extremely limited combustion within the tree pool, though set the maximum CC above zero (1%) to allow for the possibility of single trees smaller than the Landsat pixel size being consumed and not registered.

Interannual variability in KNP

By using a high-resolution FL dataset such as the one presented in this paper, it becomes possible to thoroughly investigate changes in emissions due to differing burning patterns localised to the park. To this end, we selected four individual years in the last two decades. We selected these years based on two criteria: total area burned and timing of the burns, i.e. whether the fire season was skewed towards the early or late end of the dry season. For this, we used BA data from Chen et al. (2023), available for the years 1997–2020. We excluded the years pre-2003 for consistency, as this represents the pre-Aqua MODIS era. For our comparison, we selected the years 2019 (low BA, early skew), 2013 (high BA, early skew), 2011 (high BA, late skew) and 2018 (low BA, late skew); see Fig. 4.

Fig. 4.

Burn patterns for each year from 2003 to 2020, using BA including small fires from Chen et al. (2023). All axes (including the colour bar) are relative.

GFED5 BA data were only available on a coarse grid (0.25°). As such, the ratio of burn patterns from year to year was used to adjust the Yates BA dataset for years other than 2019, rather than directly used in calculations. Emissions calculated for each year make use of variable EFs from Vernooij et al. (2023), and are presented as net emissions – that is, including only CH₄ and N₂O (Lipsett-Moore et al. 2018).

Biomass model

Results for the UAV model are shown in Fig. 5, comparing modelled with measured fuel loads for different fuel classes. R² values for surface fuels are all greater than 0.8, and for the combined total fine class, the mean percentage error is less than 20%. The model performs less well on the shrub class, R² being a somewhat lower 0.56 with a mean percentage error of just under 90%. The outputs of the S2 model are shown in Fig. 6. All classes have an R² of ~0.8. Percentage error margins are lower than in the UAV model but follow a similar pattern: lowest for total fine (9%), then litter (13%) and highest for shrubs (18%). Cumulative error for the final FL model is calculated as the sum of the root mean square errors (RMSEs) for both UAV and S2 models, and is shown in Table 4. Error is lowest for total fine fuel at ~27%, and much higher for the other classes at ~110%. Measurement error on trees in KNP is approximately 11%.

Fig. 5.

UAV RF model output compared with in situ FL measurements for wet weights of (a) grass, (b) litter, (c) total fine fuel, and (d) shrubs for both training and validation data. R² values are displayed in the top left-hand corner of each fuel class graph, along with root mean square error (RMSE) values. In parentheses are the average percentage errors. The red dotted line shows the 1:1 line. Training data are also shown as grey dots.

Fig. 6.

Sentinel-2 RF model output compared with UAV model output for wet weights of (a) total fine fuel, (b) litter (as a proxy for coarse), and (c) shrubs for both training and validation data. R² values are displayed in the top left-hand corner of each fuel class graph, along with root mean square error (RMSE) values. In parentheses are the average percentage errors. The red dotted line shows the 1:1 line. Training data are shown as grey dots.

Observational emissions

Emissions were calculated with Eqn 1 using the observational BA datasets and total DM from the UAV + S2 model. Total observed emissions for KNP in 2019 from various sources, expressed in teragrams of carbon, are shown in Fig. 7. Emissions calculated for 2019 from BA-based methods were between 5 and 6 Tg C. Pure MODIS BA-based emissions were only marginally (2%) lower than the FireCCI-based BA (4.9 and 5.1 Tg C respectively). The Yates BA-derived dataset has the highest emissions of 6.2 Tg C. Two independent emission datasets, FINN2.5 and GFAS1.2, give approximately half the carbon emissions compared with other methods for 2019 at 3.4 and 2.9 Tg C respectively. The spread of emissions is similar across all the different models – in 2019, the fire season began in March/April, peaking in June and dropping off substantially by October. The peak months (May, June, July) are where BA-based methods exceed other emissions calculation methods by as much as a factor of 2.8. BA-based emissions are largely similar in the peak months but begin to diverge later in the season (August and September). All methods show similar, low levels of biomass burning emissions for the fringe months of the fire season (March, April, October–December).

Fig. 7.

Total emissions from biomass burning in KNP from CO, CO₂, CH₄, black carbon and organic carbon for the year 2019, expressed in teragrams of carbon. Solid lines indicate emissions calculated using biomass layers generated using the UAV + S2 model, EFs from Vernooij et al. (2023) and BA products, while dotted lines indicate emissions taken directly from extant models. In parentheses in the legend is the total amount of carbon emitted in 2019 for each method.

Inter-annual variability

For coherence, in this section, we switch from reporting emissions in teragrams of carbon to teragrams CO₂ equivalent, and include four major emitted species related to biomass burning: CO₂, CO, CH₄ and N₂O. ‘Net emissions’ refers to total emissions from only CH₄ and N₂O. Of the years chosen, net emissions were highest in 2018 at 2.37 Tg CO₂ equivalent, and lowest in 2019 at 1.39 Tg CO₂ equivalent. The year with the highest BA, 2013, had marginally lower net emissions than 2018 (2.15), mostly owing to the high peak of BA in August (Fig. 8). Monthly net emissions for the year 2019 peaked in June, whereas for all other years considered, the peak fell in August. A breakdown of emissions by year and species is shown in Table 5.

Fig. 8.

Monthly emissions from KNP, showing the distribution across the year of net emissions for various burning scenarios. The total net emission, given in teragrams CO₂ equivalent, is shown in parentheses in the legend.

Fuel

The layer for which the fuel load model performs best is the mostly continuous surface fuel layer. Somewhat continuous layers with a vertical aspect (standing biomass, i.e. shrubs and trees) prove difficult to model using only 2D and height data (2D UAV or satellite-generated maps with DEMs, sometimes referred to as 2.5D data). It appears that, at least for this part of the miombo, tree or woody cover data are insufficient to get a clear picture of this biomass pool, and information about the structure of the vegetation is still necessary – such as that provided by RADAR backscattering in Bouvet et al. (2018) or via point clouds generated with LiDAR instruments (Wang et al. 2021; Qin et al. 2022). In our case, the uncertainty in the shrub class is still within a useful range, likely owing to the fact that shrubs are by definition limited in size and extent (less than 2 m), which in turn places limits on the range of biomass available in the shrub layer and creates a stronger link between surface cover and biomass – especially in savannas with broken canopy cover. This may be the main reason why it is still possible to get a reasonable AGB estimate for the shrub layer from 2.5 days’ data. With the methods applied in this paper, the tree layer proved impossible to model to a degree of accuracy that would be useful, and we conclude that proper three-dimensional structural data are necessary to capture this layer effectively. Error margins are nonetheless high for shrubs and coarse woody debris fuel classes. Although of course a lower margin of error is desirable in all cases, the coarse and shrub fuel classes are a minority contribution to total fuel consumption within KNP. Total fine fuel constitutes 60–70% of biomass burned. Coarse woody debris is the second most important fuel class, constituting less than 20% of biomass consumed. The majority of fuel consumed is therefore relatively well constrained, and the highest uncertainty in fuel load affects <25% of the total fuel consumed. In this area, future improvements could involve the use of high-resolution ground-based LiDAR to better detect small woody objects on the ground, as these proved elusive even with 5 cm resolution imagery available. This would also likely improve the shrub layer biomass with the addition of more detailed structural data. In this study, and in other studies involving larger-scale fuel consumption (e.g. van der Werf et al. 2017; van Wees et al. 2022), CC is assigned rather than measured via proxies. In our particular case, CC was assigned on the basis of a reasonable range constrained by measurements taken within KNP and moisture content. The determinants of CC in savannas include moisture content, wind speed, flame length and fuel load (Ward et al. 1996; Hély et al. 2003), though for savanna grasses at least, the direct link between moisture content and CC is strong (Hély et al. 2003). As the precipitation pattern in 2019 showed a sharp drop to practically zero from April until September, when almost all burning took place (Fig. 3), we can at least be confident that moisture content of the grasses (and likely the total fine fuel layer) for approximately 85% of the burned area was extremely low (approaching zero) and changed minimally during the season. As a result, we can be more confident that the CC parameter of this dominant fuel layer is consistently high – with some degree of uncertainty regarding the other aforementioned factors. The CC of the other layers is less certain, but the scaling approach is nonetheless one step up from larger-scale models such as GFED in the sense that we constrained CC values to ranges based on field measurements from KNP itself.

Emission models

There are three primary methods used by the fire emission models discussed in this paper to track fire activity. These are BA mapping (GFED4.1s, preliminary GFED5 and the various BA products combined with the UAV + S2 fuel layer), tracking fire activity through hotspot detection in combination with fixed FL values per biome (FINN2.5) or directly related to fire radiative energy (FRE, GFASv1.2). It is clear from Fig. 7 that the BA-based methods calculate substantially higher biomass burning emissions of total carbon for KNP. FINN2.5-derived emissions are substantially lower than GFED or UAV + S2 derived emissions as the FL in FINN2.5 is fixed across the entire southern African savanna to just over 400 g m⁻² (Wiedinmyer et al. 2023). This single value may be sufficient when considering southern African savannas as a whole, but is likely to be too low for much of KNP (see for example pink scatter dots in Fig. 5c, d, which themselves do not include tree biomass). Another possible cause of FINN2.5 showing lower emissions is the effect of using non-BA based methods for emissions calculation. Southern African savanna fires are generally faster-moving than others (Andela et al. 2019), which may limit the number of active fire detections in relation to fire size, and result in an underestimation of total burning in KNP within FINN2.5. A seemingly likely reason for the underestimation of GFASv1.2 in KNP (and also possibly elsewhere) is the tuning of GFAS to an older version of GFED (GFED3) that did not include small fires (≤100 ha). Small fires are shown to make up a significant proportion of fires worldwide, including in southern Africa (Ramo et al. 2021). However, although BA products generated on a coarser resolution (e.g. MODIS) and those on a finer grid (FireC- CISFD20) are not identical in KNP (Fig. 3), the difference is not large enough to account fully for the difference between the updated GFED versions and GFAS. One possible alternative explanation is that there is a fixed conversion factor between FRE and quantity of DM consumed for all savannas, similar to the fixed FL value in FINN. This is speculative, however, as it is hard to say how fires in KNP diverge from this value (if at all) without substantive further study. A final (and less speculative) reason for underestimation in GFAS emissions is the effect of fires beneath a canopy on measured FRE. Roberts et al. (2018) showed that a vegetation canopy can be very effective in reducing the FRE values detected at the sensor, even with broken cover. Wooded savannas have enough open space in canopy cover for visual or thermal fire detection methods (BA and active fire hotspots) to function without significant detriment even for sub-canopy fires, but FRE-based methods may suffer disproportionately in these regions. The similarity of FireCCISFD20 and MODIS-generated burned area datasets within KNP, implying a small proportion of small fires, is not unexpected. National parks or protected areas seem to be more likely to have more spatially extensive fires, owing to higher fuel connectivity and lower human footprint (Archibald et al. 2010; Eames et al. 2023). This will minimise divergence between these two BA datasets in our study area. There are few FL datasets available for the savanna region, or on a global scale, that include all fuel layers considered in this paper. An exception to this is van Wees et al. (2022), though the output of their model is DM consumed rather than a single DM layer. As this fuel consumption layer is based on MODIS BA and CC values used largely similar to those in this study, comparison of emissions generated from van Wees et al. (2022) and the UAV + S2 model combined with MODIS BA is a good (albeit indirect) method of comparing FL in both approaches. The correlation between these two datasets is extremely high and the monthly emissions are also of a very similar magnitude (Fig. 7), indicating a good agreement within the FL layer in general. However, van Wees et al. (2022) allow for a larger proportion of adult trees being consumed in fires. The fact that FL is nonetheless similar indicates that woody biomass consumption is overestimated in van Wees et al. (2022), and surface fuel underestimated within KNP. We also cannot, however, rule out at least part of this likely discrepancy coming from the fairly poorly constrained coarse woody debris layer. This incidentally also holds true for GFED4.1s and GFED5. Although both GFED versions use a slightly modified version of the MODIS BA dataset, less than 1% of the explained variance between the emissions outputs from GFED and the MODIS UAV + S2 model comes from discrepancies in the BA product used. The high level of agreement in both cases is therefore a result of a generally good match in overall FL between the two approaches, besides the likely differences in internal biomass consumption pools.

Inter-annual variability

Of the 4 years chosen, 2018 has by far the highest net emissions, a large part of which occurs in August. It seems as though the bulk of the fire season in 2018 was delayed until this point, seen in the relatively low emissions in the preceding months when compared with the other years. EFs for N₂O and especially CH₄ are higher in 2018 across the entire season (Vernooij et al. 2023), suggesting less dry and efficient burning conditions. Indeed, ERA5 data (Hersbach et al. 2020) show that not only did the 2018 wet season have the highest total precipitation (precipitation summed from July 2017 to June 2018), the bulk of this precipitation fell in the latter part of the wet season, around February. This could explain why the 2018 fire season was delayed, and also why conditions were conducive to high net emissions in that year. In contrast, precipitation in the 2018–2019 wet season was much more limited, and very little occurred in the latter portion. Conditions were likely then more conducive to burning earlier in the season, as well as more complete combustion throughout and especially in the latter part of the season. The year 2019 had the highest BA in the early dry season months (April–June) and the lowest BA overall, as well as the lowest overall net emissions. It may be that enough burning occurred in 2018 to limit the amount of BA possible in 2019; however, given the relatively low total BA for 2018, this seems less likely (see Fig. 4). Early season burning is more likely to form a patchwork of small fire scars (Ramo et al. 2021), which can prevent more extensive later season burning (Parr and Brockett 1999; Laris 2002; van Wilgen et al. 2014). Care must be taken with this conclusion from a solely emissions perspective, as applying this thinking in a wetter year such as 2018 may perhaps result in higher net emissions. The effects of savanna fire burning patterns on GHG emissions, and radiative forcing more generally, is an area where more research is needed and will be the subject of a future paper.