Coupling interactive fire with atmospheric composition and climate in the UK Earth System Model

Fire constitutes a key process in the Earth system (ES), being driven by climate as well as affecting the climate by changing atmospheric composition and impacting the terrestrial carbon cycle. However, studies on the effects of fires on atmospheric composition, radiative forcing and climate have been limited to date, as the current generation of ES models (ESMs) does not include fully atmosphere–composition–vegetation coupled fires feedbacks. The aim of this work is to develop and evaluate a fully coupled fire–composition–climate ES model. For this, the INteractive Fires and Emissions algoRithm for Natural envirOnments (INFERNO) fire model is coupled to the atmosphere-only configuration of the UK's Earth System Model (UKESM1). This fire–atmosphere interaction through atmospheric chemistry and aerosols allows for fire emissions to influence radiation, clouds and generally weather, which can consequently influence the meteorological drivers of fire. Additionally, INFERNO is updated based on recent developments in the literature to improve the representation of human and/or economic factors in the anthropogenic ignition and suppression of fire. This work presents an assessment of the effects of interactive fire coupling on atmospheric composition and climate compared to the standard UKESM1 configuration that uses prescribed fire emissions. Results show a similar performance when using the fire–atmosphere coupling (the “online” version of the model) when compared to the offline UKESM1 that uses prescribed fire. The model can reproduce observed present-day global fire emissions of carbon monoxide (CO) and aerosols, despite underestimating the global average burnt area. However, at a regional scale, there is an overestimation of fire emissions over Africa due to the misrepresentation of the underlying vegetation types and an underestimation over equatorial Asia due to a lack of representation of peat fires. Despite this, comparing model results with observations of CO column mixing ratio and aerosol optical depth (AOD) show that the fire–atmosphere coupled configuration has a similar performance when compared to UKESM1. In fact, including the interactive biomass burning emissions improves the interannual CO atmospheric column variability and consequently its seasonality over the main biomass burning regions – Africa and South America. Similarly, for aerosols, the AOD results broadly agree with the Moderate Resolution Imaging Spectroradiometer (MODIS) and the Aerosol Robotic Network (AERONET) observations.

1 Introduction

Fires can exert a substantial forcing on the Earth's climate by affecting different components of the Earth system (ES) such as the biosphere, atmosphere and cryosphere (Bowman et al., 2009; Daniau et al., 2013). Changes in vegetation cover caused by the fire modifies the regional to local-scale surface albedo, soil water holding capacity and surface evaporation, resulting in complex interactions and feedbacks within the climate system (Li et al., 2017; Myhre, 2005). In addition, fire emissions contribute to the global budgets of greenhouse gases (methane, ozone) and aerosol particles (black carbon, organic carbon) (Lasslop et al., 2019; Voulgarakis and Field, 2015), resulting in direct and indirect effects on solar irradiation as well as changes in the land surface by means of black carbon deposition, which, in turn, leads to modifications of the surface albedo of bright ice and snow surfaces (Ramanathan and Carmichael, 2008; Thomas et al., 2017). Moreover, climate variability and climate change can also impact fire frequency and other aspects of fire behaviour. Among others, Gillett et al. (2004) and Westerling et al. (2006) presented evidence that climate change has contributed to an increase in fire frequency in North America and Eurasia. However, a long-term increase in the length of the fire season or in weather conditions conducive to wildfires does not necessarily lead to an increase in burned area, as this is also limited by the available fuel (Doerr and Santín, 2016). These previous studies rely on statistical models of fire danger and burned area forced with several different climate projections and do not account for the composition–climate feedback effects and interactions caused by fire emissions, changes in vegetation productivity and structure or fire–vegetation–climate interactions (Rabin et al., 2017). Thus, the effort towards being able to represent fire within the ES framework is of great importance, particularly since fire as a process is highly coupled within the ES and responds to both natural and anthropogenic changes. Comprehensive model representations of the ES should therefore include the effects of fires on the climate and the effects of climate on fires (Rabin et al., 2017). In addition, it is important to consider not only past and present climate and different future climate scenarios but also scenarios of demographic changes since human population size, distribution and economic activity play a role in the occurrence of fires and amount of burnt area but have so far received less attention.

Recent studies have contributed to the understanding and quantification of various aspects of the effects of fires on climate. Some of these studies have focused on specific impacts and/or fire events (e.g. López-Saldaña et al., 2015; Samset et al., 2014; Baker et al., 2016; Bali et al., 2017; Dintwe et al., 2017; Li et al., 2017), while others have looked at fire emissions (e.g. Voulgarakis and Field, 2015; Baker et al., 2016; Lasslop et al., 2019) or fires that occur within a particular ecosystem or region (e.g. Hirota et al., 2011; Athanasopoulou et al., 2014; Rogers et al., 2015; Dintwe et al., 2017). On the other hand, global-scale assessments highlight the complexity and uncertainties of these impacts, particularly those from aerosols, as well as the difficulty in performing a comprehensive analysis at a global scale accounting for all relevant levels of interactions (Huang et al., 2016; Unger and Yue, 2014). As such, the total radiative effect of fires remains fundamentally uncertain, making climate–fire feedbacks relevant in the context of climate change research (Carslaw et al., 2010; Unger and Yue, 2014; Ward et al., 2012).

Despite the importance of climate–fire feedbacks and the large impacts of fire on the ES and its net radiative forcing, there is still a large knowledge gap in this subject. As has been pointed out by the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (Settele et al., 2015), this, together with the varying projections of future climate leads to large uncertainties in the direction of regional changes in future fire regimes. Still, the work from Ward et al. (2012) provides a notable global overview of the most important effects while maintaining a consistent method of analysis across the ES by providing a model-based analysis that examines how the radiative forcing has changed since pre-industrial times while also providing an estimate on how much it may change in the future. In their study, the authors compared a set of experiments with fire emissions with a no-fire emissions scenario and found that fires have an overall net negative global radiative forcing of approximately −1.02 W m⁻² across the study periods, centred around 1850 to 2000 and 1850 to 2100. The authors also found that the masking of fire aerosol impacts on clouds by anthropogenic aerosols between 1850 and 2000 decreases the magnitude by 0.6 W m⁻² of the fire radiative forcing for that historical period. According to Kloster et al. (2012) and Pechony and Shindell (2010), for the 2000–2100 period, global emissions from fires primarily depend on the climate forcing. Ward et al. (2012) has also shown that even though models may have an overall similar net change in radiative forcing over time, this net forcing can be driven by different sources – background anthropogenic and natural (e.g. dust, sea salt) emissions or biomass burning emissions – which balance each other to obtain the same net result. Therefore, the choice of model also proves to be a source of uncertainty in the existing published results. Voulgarakis et al. (2015) have shown that fires play a large and even dominant role in driving the interannual variability of key trace gases and aerosols that influence air quality and climate. For example, fires are almost entirely responsible for the interannual variability of carbon monoxide and carbonaceous aerosols, with a major role also in hydroxyl radical (OH) interannual variability.

The large uncertainties in net radiative forcing of fires are predominantly caused by the uncertainties in the total fire emissions, their spatial variability, model representation of aerosol–cloud interactions, the simulated land cover and future atmospheric composition. Fires are also the largest source of carbonaceous aerosol globally accounting for 55 % to 60 % of the of primary organic carbon (OC) and black carbon (BC) aerosol emissions and are the dominant source of aerosol emissions for the central African and the Amazonian basin regions (Andreae and Rosenfeld, 2008; Mahowald et al., 2011; Ward et al., 2012). Furthermore, there remain knowledge gaps and a need for improvement in model estimates of the distribution of fires (Kloster et al., 2012). Some progress has already been made in recent years by developing global fire emission inventories, which are essential for both model development and validation (Vongruang et al., 2017; Van Der Werf et al., 2017). Equally crucial is having a robust representation of fires within land surface models which has been found to be inaccurate to varying degrees. For example, very few fire models include the representation of peatland fires, which have been estimated to represent an average emission source of ∼100 to 200 Tg C yr⁻¹, for the period between (1997–2009), 10 % of the global total carbon fire emissions (∼2.0 Pg C yr⁻¹) (Van Der Werf et al., 2010), and dominate global fire emissions variability (Voulgarakis et al., 2015; Van Der Werf et al., 2010). Over the last decade, fire modelling has seen important advances led mainly by the interest of the research community to incorporate these into Earth system models (ESMs) with the aim of studying fire–climate–composition interactions in a fully coupled fashion.

The goal of this work is to develop and evaluate a fully coupled fire–climate–composition ES model. For this, we have built on the work of Mangeon et al. (2016) on the INteractive Fires and Emissions algoRithm for Natural envirOnments (INFERNO) and have coupled this fire model to the atmospheric component of the UK Earth System Model, version 1 (UKESM1; Sellar et al., 2019). This new fire–atmosphere interaction replaces the prescribed transient monthly varying biomass burning emissions in UKESM1. As a result, through atmospheric chemistry and aerosols, the interactive fire emissions can affect radiation and clouds, thereby affecting weather/climate and the meteorological drives of fires themselves.

The INFERNO model, as described by Mangeon et al. (2016), excludes any representation of socio-economic factors that influence both fire ignition and suppression. Several studies have used the gross domestic product (GDP) as a proxy to represent the role of socio-economic factors influencing fires (Aldersley et al., 2011; Bistinas et al., 2014; Hantson et al., 2016). However, GDP does not account for socio-economic policies and factors that can also impact the management of fire (Ganteaume et al., 2013; Pausas and Keeley, 2014; Pezzatti et al., 2013). With this in mind, INFERNO was updated to improve the representation of population dynamics and economic activity for anthropogenic ignition and suppression of fire. This work presents an assessment and evaluation of the effects of interactive fires on atmospheric composition and climate compared to the standard UKESM1 configuration, which so far has used prescribed fire emissions.

2 Methods and data

2.1 Fire model – INFERNO

The INFERNO model developed by Mangeon et al. (2016) is the integrated fire model for the Joint UK Land Environment Simulator (JULES; Best et al., 2011; Clark et al., 2011), which serves as the land surface component of UKESM1. INFERNO uses an approach based on the work of Pechony and Shindell (2009) adapted to allow interaction within an ESM framework. More precisely, water vapour pressure deficit is used as one of the main indicators of biomass flammability in the model, while an inverse exponential relationship is used to relate flammability to soil moisture. In INFERNO, fire ignitions can be caused by cloud-to-ground lightning strikes and anthropogenic ignitions in the following three ways: derived from a multi-year annual mean; assuming constant human ignitions but with varying cloud-to-ground lightning strikes, which always strike to start a fire and therefore accounts for natural variability in fire ignitions, which can be simulated interactively with an ESM or prescribed from observations; using varying human ignitions together with natural ignitions as described in the second mode.

The burnt area for a plant functional type (PFT), (BA_PFT) (fraction s⁻¹), is given by Eq. (1):

$$\begin{array}{}\text{(1)}& {\text{BA}}_{\text{PFT}}=I_{\mathrm{T}}{F}_{\text{PFT}}\stackrel{\mathrm{‾}}{{\text{BA}}_{\text{PFT}}},\end{array}$$

where I_T represents the fire ignitions ($\mathrm{ignitions}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), including natural and human ignitions as well as fire suppression by humans, F_PFT is the flammability per PFT dependent on the 1.5 m temperature, 1.5 m relative humidity and fuel density as defined in Eq. (4) through Eq. (6) from Mangeon et al. (2016), and $\stackrel{\mathrm{‾}}{{\text{BA}}_{\text{PFT}}}$ (m² ignition⁻¹) is the scaled average burnt area per ignition for each PFT, which decouples the fire spread stage from local meteorology and topography, processes which are not typically resolved in coarse grids, such as the those often used in ESMs. The values used here are an adaptation of those reported by Andela et al. (2018) to the model PFT setup as shown in Table 1.

Table 1 Biomass burning scaled average burnt areas (km² fire⁻¹) and emission factors (g kg⁻¹ dry matter) for species emitted from various plant functional types based on Andela et al. (2018) and Andreae (2019).

The emitted carbon per PFT, EC_PFT ($\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), in Eq. (2) is calculated based on the burnt area and combustion completeness, accounting for the wetness of fuel. INFERNO also represents emissions for trace gases: carbon dioxide (CO₂), carbon monoxide (CO), methane (CH₄), nitrogen oxide (NO_x), sulfur dioxide (SO₂), ethane (C₂H₆), propane (C₃H₈), formaldehyde (HCHO), Acetaldehyde (MeCHO), acetone (Me2CO), ammonia (NH₃), dimethyl sulfide (DMS); and aerosols: OC and BC. These emissions are estimated based on INFERNO's emitted carbon estimate (Eq. 2) by using the emission factors based on the work by Andreae (2019) as shown in Table 1.

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\text{EC}}_{\text{PFT}}=& {\text{BA}}_{\text{PFT}}\sum _{\text{leaf},\text{stem}}^i\\ & \left[{\text{CC}}_{\text{min},i}+\left({\text{CC}}_{\text{max},i}-{\text{CC}}_{\text{min},i}\right)\left(\mathrm{1}-\mathit{\theta }\right)\right]{\mathrm{C}}_i\end{array}\end{array}$$

CC_min and CC_max are the minimum and maximum combustion completeness for both leaves (CC_min=0.8 and CC_max=0.9) and stems (CC_min=0.2 and CC_max=0.4), C_i is the carbon stored in each PFT leaves or stem and θ the unfrozen soil moisture as a fraction of saturation.

INFERNO fire ignitions are split into natural ignitions (I_N) (ignition m² s⁻¹) from cloud-to-ground lightning (externally provided to INFERNO) and anthropogenic ignitions (I_A) (ignition m² s⁻¹) which are dependent on population density (PD) (people m⁻²) as described in Eq. (3). Moreover, it is assumed that humans are also responsible for suppressing fires which is accounted for through a suppression function dependent on human population density given by Eq. (4). Originally, Eqs. (3) and (4) were developed to include only information on population density. We now include a human development index (HDI) term (1−HDI) in these equations that represents socio-economic factors impacting fire ignition and suppression. HDI is calculated based on three indicators designed to capture the income, health and education dimensions of human development. For areas where there is more effort in human development improvements, it is assumed that fire ignition is decreased, and fire suppression is increased.

HDI data were obtained from the gridded global datasets for gross domestic product and human development index (Kummu et al., 2018). The anthropogenic ignitions (I_A) (ignition m² s⁻¹) is represented by Eq. (3), the fraction of fires not suppressed by humans (f_NS) by Eq. (4) and the total ignitions (I_T) are represented by Eq. (5).

$$\begin{array}{}\text{(3)}& {\displaystyle }{I}_{\mathrm{A}}=& {\displaystyle }{k}_{\text{(PD)}}\text{PD}\mathit{\alpha }\times \left(\mathrm{1}-\text{HDI}\right)\text{(4)}& {\displaystyle }{f}_{\text{NS}}=& {\displaystyle }\mathrm{7.7}\left(\mathrm{0.05}+\mathrm{0.9}\times e^{-\mathrm{0.05}\text{PD}}\right)\times \left(\mathrm{1}-\text{HDI}\right)\text{(5)}& {\displaystyle }{I}_{\mathrm{T}}=& {\displaystyle }\left(I_{\mathrm{N}}+I_{\mathrm{A}}\right){\displaystyle \frac{f_{\text{NS}}}{\mathrm{8.64}\times {\mathrm{10}}^{\mathrm{10}}}},\end{array}$$

where $k_{\text{(PD)}}=\mathrm{6.8}\times \text{PD}-\mathrm{0.6}$ is a function that represents the varying anthropogenic influence on ignitions in rural versus urban environments, the parameter α=0.03 represents the number of potential ignition sources per person per month per km², I_N represents the natural ignitions due to lightning and HDI the human development index.

Furthermore, it should be highlighted that in this configuration of INFERNO, there are no interactions between fire and vegetation, and it does not include a peat-burning capability.

2.2 UK Earth System Model (UKESM1)

This study uses the atmospheric and land components of UKESM1 (Sellar et al., 2019, 2020) following the protocol set by the Atmospheric Model Intercomparison Project (AMIP, Eyring et al., 2016). The model resolution used in this configuration is N96L85. This is equivalent to a horizontal resolution of 135 km in the midlatitudes and 85 terrain-following vertical levels ranging up to an altitude of 85 km above sea level. The science configuration of the atmosphere component is based on the Global Atmosphere 7.1 (GA7.1) and the Global Land 7.0 (GL7.0) as described by Walters et al. (2019) used in the configuration of the Hadley Centre Global Environment Model version 3 (HadGEM3; Hewitt et al., 2011) coupled to the terrestrial carbon/nitrogen cycles (Sellar et al., 2019) and interactive stratosphere–troposphere chemistry (Archibald et al., 2020) from the UK Chemistry and Aerosol (UKCA; Morgenstern et al., 2009; O'Connor et al., 2014) model. The framework also includes the UKCA prognostic aerosol GLOMAP-mode scheme (Carslaw et al., 2010; Mulcahy et al., 2019), where secondary aerosol formation is determined by interactive oxidants from the UKCA stratosphere–troposphere chemistry scheme (Archibald et al., 2020).

As per the AMIP protocol, sea surface temperature and sea ice are taken from the unmodified dataset of Durack and Taylor (2017) and horizontally interpolated to the model resolution. In this model setup, the dynamic vegetation model (TRIFFID, Cox, 2001) is deactivated and replaced by prescribed vegetation properties from a coupled historical simulation with the same base model, as shown in Fig. 1, preserving consistency in the forcing due to land use change between the UKESM1 coupled and AMIP experiments. In a similar fashion, seawater concentrations of dimethyl sulfide (DMS) and chlorophyll-a monthly climatologies are taken from the coupled historical experiment and are used by the atmosphere model top calculates fluxes of DMS and primary marine organic aerosol (Mulcahy et al., 2019).

Figure 1 Dominant vegetation PFT per grid box for as prescribed in UKESM1-AMIP configuration.

External forcing datasets for biomass burning aerosol and trace gas emissions follow those stipulated under the coordination of the CMIP6 protocol (Van Marle et al., 2017). These are a combination of satellite observations from 1997 with various proxies and modelling results from six fire models which have participated in the Fire Model Intercomparison Project (Rabin et al., 2017) aimed at providing a dataset of biomass burning emissions for use in CMIP6. UKESM1 uses emissions of primary carbonaceous aerosol – black carbon (BC) and organic carbon (OC) – as well as DMS, which acts as a precursor to sulfate aerosol. In terms of gas-phase biomass burning emissions, UKESM1 uses lumped emissions of ethene (C₂H₄) and ethane (C₂H₆) as C₂H₆, and emissions of propane (C₃H₈), formaldehyde (HCHO), acetone (CH₃)₂CO), acetaldehyde (CH₃CHO), carbon monoxide (CO) and nitric oxide (NO). Aerosol emissions of BC and OC from biomass burning are spread evenly in the vertical over the first 20 model levels – corresponding to the lowest 3 km – and are treated with a geometric mean diameter of 150 nm, while all other biomass burning trace gas species are injected into the model's lowest layer and mixed simultaneously by the boundary layer mixing scheme. A biomass burning emissions scaling factor of 2 is applied to biomass burning aerosols, following the evaluation work of Johnson et al. (2016), which improves the agreement between observed and simulated aerosol optical depth (AOD) across the three evaluated wave lengths (440, 550, 700 nm) when compared to observations. This AOD bias is also evident in the comparison of previously published top-down and bottom-up estimates and using a scaling factor correction is a practice widely used in the modelling community due to current model deficiencies (Kaiser et al., 2012).

Annual anthropogenic emissions of reactive gases are prescribed to the model and are taken from the Community Emissions Data System (CEDS, Hoesly et al., 2018) as prepared for use in CMIP6. Both biomass burning and anthropogenic emissions datasets are regridded from their native resolution to N96L85 while conserving global annual totals and seasonal cycles.

A set of natural emissions of other species which are not simulated by UKESM1 is prescribed through precomputed fluxes. This includes emissions of oceanic emissions of CO, C₂H₆ (including C₂H₄ lumping) and C₃H₈ (including propene (C₃H₆) lumping) from POET (Granier et al., 2005) and correspond to the annual cycle inventory for the year 1990 (12 monthly fluxes). Biogenic emissions of CO, HCHO, MeOH, C₃H₆ and C₃H₈, as well as CH₃CHO and MeCHO are taken from the MACCity-MEGAN emissions inventory (Sindelarova et al., 2014) and are provided to the model based on the 2001–2010 monthly mean climatology. Soil emissions of NO_x are distributed according to Yienger and Levy (1995) and scaled to give a global annual total of 12.0 Tg NO yr⁻¹, perpetually applied to all years.

Further details on the UKESM1 model setup can be found in Sellar et al. (2020) and further details (e.g. lightning NO_x emissions, lower boundary conditions, interactive BVOC emissions) and an evaluation of the performance of the UKCA chemistry and aerosol schemes in UKESM1 are available in Archibald et al. (2020) and Mulcahy et al. (2019).

2.3 Fire–composition–atmosphere coupling

JULES is the land surface model used in UKESM1, and a coupling interface is in place for the exchange of variables and drivers between the atmosphere and land components. For this work, the interface was extended to allow the coupling of the required atmospheric variables from the atmospheric model to INFERNO through JULES. The atmospheric model provides to INFERNO the surface pressure, 1.5 m temperature, 1.5 m specific humidity and precipitation at every model time step. Temporal variations of lightning can have a large impact on the simulation of burned area (Felsberg et al., 2018). With this in mind, and to maintain consistency between fire lighting ignitions, the state of the atmosphere and atmospheric composition, the cloud-to-ground lightning simulated by UKESM1 is passed down to INFERNO. The lighting parametrization used follows the work of Price and Rind (1994), which makes use of parameterized lightning flash frequency of $\mathrm{3.44}\times {\mathrm{10}}^{-\mathrm{5}}{H}^{\mathrm{4.9}}$ min⁻¹ over land and $\mathrm{6.4}\times {\mathrm{10}}^{-\mathrm{4}}{H}^{\mathrm{1.73}}$ min⁻¹ over ocean (where H is the cloud depth in kilometres), along with a cloud–cloud and cloud–ground flash ratio based on the grid-cell latitude. The cloud depth is determined using the convective cloud base and top levels diagnostics from the convection scheme.

Conversely, in order to allow for fire–climate interactions through atmospheric chemistry and aerosols, the coupling framework was extended to pass INFERNO-derived emissions to the UKCA chemistry and aerosol model of CO, NO_x, C₂H₆, C₃H₈, HCHO, MeCHO, Me2CO, NH₃, DMS, OC and BC at every model time step. As described by Archibald et al. (2020), in UKESM1, biomass burning emissions of C₂H₆ and C₃H₈ include emissions of C₂H₄ and C₃H₆, respectively. In this interactive fire emissions framework, we do not consider this emission aggregation, and considering the small contribution of biomass burning to these species (15.7 % for C₂H₆ and 13.6 % for C₃H₈), we do not consider this to have a large impact on the modelled results.

Aerosol emissions are distributed vertically following an exponential increasing function (emission increasing with height) from the first model level to a fixed top of plume height defined as the 20th model level (∼3 km). This approach is a simplification based on the work of Rémy et al. (2017). This new fire–atmosphere interaction replaces the prescribed transient monthly varying biomass burning emissions in UKESM1. As a result, through atmospheric chemistry and aerosols, the interactive fire emissions can affect radiation and clouds, thereby affecting weather/climate and the meteorological drives of fires themselves.

The model experiments were run for the period 1974 to 2014 and the initial 5 years (1974–1979) were considered as simulation spin-up and discard from the analysis.

For the sake of brevity, from here onwards, this fire–composition–atmosphere coupled configuration based on UKESM1 is referred to as UKESM1+INFERNO.

2.4 Burnt area and emissions evaluation

Burnt area data from the Global Fire Emissions Database version 4 (GFED4s) (Giglio et al., 2013) were used to assess the performance of the UKESM1+INFERNO in modelling fires. This dataset is provided as a gridded product at a 0.25^∘ resolution and is derived from a multi-sensor satellite dataset, including satellite dataset based on active fire detection, including small fires, based on statistical modelling, as detailed in Randerson et al. (2012).

Figure 2 Basis regions, as defined in the GFED4s dataset (Giglio et al., 2013).

The basis regions as defined in the GFED4s dataset (Fig. 2) were applied to the modelled data as required in order to perform a regional assessment of the INFERNO results. From here after, these regions will be named according to the acronyms defined in Fig. 2.

To evaluate the fire emissions from INFERNO within the global ESM, the data from Global Fire Assimilation System (GFAS) were used. GFAS calculates biomass burning emissions by assimilating satellite observations of fire radiative power (FRP) from the Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua and Terra satellites (Kaiser et al., 2012). Combustion rates in GFAS are calculated with land-cover-specific conversion factors which have been derived from a linear regression analysis between the FRP of GFAS and the dry matter combustion rate of GFED (Heil et al., 2010). Emissions for 40 gas-phase and aerosol trace species are then calculated by applying emission factors based on the works of Andreae and Merlet (2001), Christian et al. (2003) and Akagi et al. (2011). Daily emissions are available on a global $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid from 2003 to the present day.

2.5 AOD and extinction coefficient evaluation

Both ground-based measurements and satellite retrievals of AOD are used to evaluate the model performance. Ground-based measurements from the Aerosol Robotic Network (AERONET) provides quality assured measurements of aerosol optical properties across the globe (Holben et al., 1998, 2001). For this study, the monthly AOD at 440 nm is used from the version 2 level 2.0 product for stations that have at least a 5-year overlapping period with the model simulations. AOD evaluation is complemented with data from satellite retrievals. The MODIS aerosol product provides daily observations of the AOD with a global coverage. We used the collection of the level 3 V6 MODIS monthly data of AOD at 550 nm for the period from January 2003 to December 2012, which is produced from daily means blending the Dark Target and Deep Blue algorithms (Hsu et al., 2004; Sayer et al., 2014) and provided as gridded data at a 1^∘ spatial resolution.

To evaluate the vertical profile of aerosol, the Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) lidar level 3 aerosol extinction coefficient profiles product at 532 nm were used. This product reports monthly mean profiles of aerosol optical properties and is quality screened prior to averaging on a uniform spatial grid with a meridional and zonal resolution of 2 and 5^∘ respectively and a vertical resolution of 60 m. For comparison with model results, the focus is on clear sky averages (Tackett et al., 2018).

2.6 Carbon monoxide

Evaluation of the gas-phase biomass burning emissions is focused on carbon monoxide (CO), the most abundant chemically active pollutant emitted by fires (van der Werf et al., 2010). For this, model results are evaluated against observations from the Tropospheric Emission Spectrometer (TES-AURA), which have been used before to evaluate CO in previous versions of UKCA (O'Connor et al., 2014; Voulgarakis et al., 2011). We have used TES-AURA-Lite version 007 data (Beer et al., 2001; Bowman et al., 2006) covering the 5-year period between 2007–2012 where there were fewer data gaps. To compare model results with TES-AURA observations, the hourly CO model output has been interpolated onto the 14 TES-AURA-Lite pressure levels, as well as satellite swath location. Furthermore, the TES-AURA sampling and averaging kernels and a priori profiles were also applied to the model data. Both TES-AURA and modelled processed data are then monthly averaged into a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ resolution grid. Only the vertical region where TES-AURA is more sensitive to CO was used – average over the vertical column between 700 and 300 hPa.

3 Results

In this section, the results of the implementation of the fire–composition–atmosphere coupling in UKESM1+INFERNO are analysed. When comparing datasets (model or observed) with different grid resolutions, the higher resolution dataset is regridded to the lowest resolution grid using a first-order conservative area-weighted regridding method. Statistical significance of the differences presented here were examined using a Student's t test (Wilks, 2011) with a 95 % confidence level.

3.1 Burnt area

Figure 3 shows the annual mean burnt area fraction (1997–2010) for (a) UKESM1+INFERNO and (b) GFEDv4. The overall geographical pattern of the annual average burnt area fraction is well reproduced by the model with a global pattern correlation of 55.3 % when compared with GFEDv4. The model represents the observed pattern in the major fire regions: South America, Africa and Eurasia. On the other hand, the model underestimates the northern Australia fires and boreal regions. Although the burnt area fraction over Africa is well represented, there is a large (50 %) underestimation of the fires over Africa in the northern African region (NHAF). This underestimation can be attributed to the Saharan bare soil extending too far south, causing a lack of grassland in the Sahel region, which is a result of precipitation deficits associated with errors in the position and intensity of monsoon systems (Sellar et al., 2019; Williams et al., 2018). In addition, there is an overestimation of tree fraction in savanna biomes, such as the southern African region (SHAF) and the southern edge of the Amazon rain forest region (SHSA). The differences in the specified scaled average burnt areas for these biomes – smaller for trees than grasses – cause an underestimation of fire size in these regions. This overestimation of the tree fraction is attributed to the lack of fire disturbance in the fully coupled UKESM1 configuration, the inclusion of which could potentially improve vegetation structure in these regions (Burton et al., 2019).

Figure 3 Burnt area fraction (%) mean annual average (1997–2010) for (a) GFED4s and (b) UKESM1+INFERNO.

Figure 4 Burnt area annual mean time series (a, b) anomaly relative to climatology (Mha) and (c) burnt area fraction climatology (%) for UKESM1+INFERNO (solid line) and GFEDv4 (dashed line) with the shaded area representing the standard error around the GFEDv4 data.

These biases found for the different regions have an impact on the modelled global fire behaviour. Although there is a global underestimation of the annual average burnt area of approximately 250 Mha (Fig. 4a), the model is able to capture the negative trend found in the observations (Fig. 4b). In terms of seasonality, the model produces a bimodal pattern in the fire activity, as observed in GFEDv4, which peaks around both the austral and boreal late summer season.

Although there is a global underestimation contributing to the interannual variability of burnt area, it is important to stress that the biases found in NHAF and SHAF contribute to most of the biases found in the global burnt area climatology. A November-to-January difference of 45 Mha for the NHAF region combined with a July-to-September difference of 30 Mha results in an absolute bias of greater than 300 %, as can be seen in Fig. 5.

Figure 5 Burnt area fraction climatology (%) for UKESM1+INFERNO (solid line) and GFEDv4 (dashed line) with the shaded area representing the standard error around the GFEDv4 data for NHAF (a) and SHAF (b).

3.2 Biomass burning emissions

In order to develop a full fire–composition–climate ES model, it is paramount to evaluate the ability of the INFERNO – UKESM1 coupling with regards to providing realistic emissions of chemistry and aerosol species. Currently this coupling framework provides trace gases emissions of CO, NO_x, C₂H₆, C₃H₈, HCHO, DMS, NH₃ and DMS, as well as aerosol emissions of OC and BC. These emissions are estimated based on INFERNO emitted carbon by using PFT-specific emission factors. For this reason, the different species present similar broad characteristics. With this in mind and considering the observed datasets available and used to assess the model performance, we will focus on the analysis of CO, OC and BC.

Figure 6 Biomass burning emissions mean annual average (2003–2011) for CO (g m⁻²) in (a) GFAS, (d) UKESM1+INFERNO and (g) difference between UKESM1+INFERNO and GFEDv4; OC (g m⁻²) in (b) GFAS, (e) UKESM1+INFERNO and (h) difference between UKESM1+INFERNO and GFEDv4 and BC (g m⁻²) in (c) GFAS, (f) UKESM1+INFERNO and (i) difference between UKESM1+INFERNO and GFAS. Differences are only shown for statistically significant points with a 95 % confidence level.

As seen for the burnt area results, the overall global pattern for annual average biomass burning emissions is, in general, well reproduced by the model (Fig. 6). The main regions of biomass burning emissions, Africa and South America, are captured in the global spatial pattern. However, there is a large overestimation of the biomass burning emissions, despite underestimation of area burned, in the southern edge of NHAF, SHAF as well as the eastern side of SHSA for all the species (difference>300 %). Furthermore, in the SHAF region, the emissions extend further south into the midlatitudes. These overestimations can be attributed to the overestimation of the tree fraction in these regions leading to a greater content of carbon available for combustion. On the other hand, the model underestimates the emissions for the boreal regions (emissions close to zero in the model). These biases result in smaller values of global pattern correlation for biomass burning emitted species of 36.7 % for CO, 38.4 % for OC and 53.2 % for BC when compared with GFAS. UKESM1+INFERNO shows a good climatological performance (no large bias; Fig. 6) for equatorial Asia (EQAS). However, it fails to reproduce specific large fire events (that are infrequent) associated with peatland fires that represent a substantial amount of the global biomass burning emissions (not shown). We do not expect the model to be able to capture such emission events, as the current land surface model does not include a peatland land class, and therefore we do not include a peat-burning capability. However, it should be noted that although we have a negative bias relative to the observations, the observations themselves may be biased due to the lacking efficiency in detecting low-temperature smouldering peat fires in MODIS products used by GFED4s and GFAS.

Figure 7 Annual mean time series and seasonal cycle climatology of biomass burning emissions (mg m⁻²) of CO (a, b), OC (c, d) and BC (e, f) for UKESM1+INFERNO (solid line) and GFAS (dashed line). The shaded area represents the standard error in the GFAS data.

The biomass burning emission (kg m⁻²) annual mean time series and climatology (Fig. 7) show that there are large biases in the annual mean time series for all the emissions. When compared to GFAS, CO shows a root mean squared error (RMSE) of 349.53 mg m⁻² and a correlation of 74.79 %, whereas OC shows a RMSE of 10.08 mg m⁻² and BC a RMSE of 1.08 mg m⁻². Both simulated OC and BC show non-significant correlation (when tested with a 95 % confidence level) of 56.06 % and 49.41 %, respectively. Despite this, there is a good agreement to their observations with regards to the interannual variability and climatology, with the model capturing the negative emission trends in the 1997–2010 period and being able to reproduce the observed seasonal cycle – higher emissions of CO, OC and BC during the period from June to October (in Fig. 7b, d and f, respectively). Despite the non-significant correlation between the interannual time series, there is good qualitative agreement for both interannual and seasonal variability. This is partly due to the compensation of regional biases and partly due to the burnt area bias. The low correlations obtained for OC and BC interannual time series are mainly due to the overestimation of the negative trend by the model.

3.2.1 Sensitivity to land surface

In order to better understand the impact of the overestimation of tree fraction in savanna biomes, sensitivity experiments were performed where the dominant tree PFT is replaced with an equal quantity of the dominant grass PFT. We focus on a region over SHAF, characterized by a savanna biome, but where in UKESM1 is dominated by broadleaf evergreen tropical trees (∼70 % of total PFT fraction) – latitudes between 3.0 and 29.0^∘ S and longitudes between 7.0 and 42.0^∘ E.

Two sensitivity experiments were performed where, for this specific region, the broadleaf evergreen tropical trees are replaced by C₄ grasses and run for the period 1980–1985:

• T2G10 – 10 % of trees are changed to grass.

• T2G50 – 50 % of trees are changed to grass (this experiment has the closest representation of the observed land surface).

In the SHAF region, which includes the area where the PFTs were changed; there are 12.17 % and 79.83 % increases of the burnt area for 10 % and 50 % changes in the tree cover into grasses, respectively (T2G10 and T2G50, Table 2). Thus, regionally, the modelled burnt area is hypersensitive to a change of the underlying vegetation (the change in burnt area is greater than the change in vegetation cover fraction). In contrast, for the case of biomass burning emissions, there is much less sensitivity to a change in the vegetation cover. In this case, for a 10 % change of the dominant tree to the dominant grass PFT (T2G10), there were no statistically significant changes to emissions. However, a statistically significant decrease of ∼14 %–25 % is found for T2G50, showing that, although biomass burning emissions are hyposensitive to a change of the underlying vegetation, when there is a larger change to the underlying vegetation, this can result in significant changes.

Table 2 Total annual average values and relative change (%) compared to UKESM1+INFERNO (period between 1980–1985) of burnt area (Mha), carbon monoxide (g m⁻²), organic carbon (g m⁻²) and for black carbon (g m⁻²) for UKESM1+INFERNO and the different sensitivity experiments, T2G10 and T2G50 – global and SHAF regions. Values in bold show that the difference between a given experiment and UKESM1+INFERNO is statistically significant with a 95 % confidence level.

This sensitivity to the underlying vegetation can also have a significant impact at the global scale for both burnt area (statistically significant increase of 5.05 % and 19.00 % for T2G10 and T2G50, respectively) and biomass burning emissions (significant decreases of 5.26 % and 8.33 % for OC and BC in T2G50, respectively). The change applied to the vegetation causes significant changes in both burnt area and biomass burning emissions not only locally but also for regions away from the area where the vegetation was changed, through changes in atmospheric dynamics at a global scale (Fig. 8). However, further work is required to assert if these changes are mainly associated with the change to the land cover or the fire–atmosphere–composition feedbacks. This strongly indicates that a realistic representation of the vegetation distribution could significantly improve the model performance when compared to observations.

Figure 8 Total annual average fractional of change (period between 1980–1985) of burnt area (Mha), carbon monoxide, organic carbon and for black carbon for the UKESM1+INFERNO and the different sensitivity experiments, T2G10 and T2G50. Stippling is shown for points where the difference between a given experiment and UKESM1+INFERNO is statistically significant with a 95 % confidence level. The green rectangle indicates the area where the prescribed vegetation was changed.

3.3 Carbon monoxide atmospheric column

When analysing the CO mean column mixing ratio averaged between 700 and 300 hPa (Fig. 9), it is possible to see that the atmospheric column of CO is dominated by the hot spots of biomass burning over South America and Africa and anthropogenic emissions in the Northern Hemisphere, with a strong north–south hemispheric gradient due to the short life time of CO – approximately 30 d – compared to the timescale of interhemispheric mixing.

Figure 9 Column volume mixing ratio carbon monoxide (ppm) averaged between 700 and 300 hPa over the period 2007–2012 for (a) TES-AURA satellite retrievals, (b) UKESM1+INFERNO, (c) difference between UKESM1 and TES-AURA, and (d) difference between UKESM1+INFERNO and TES-AURA.

When comparing the model results to TES-AURA, it can be seen that there is an underestimation of the column CO in the northern hemisphere. This underestimation happens both in UKESM1 (Fig. 9c) and UKESM1+INFERNO (Fig. 9d). As documented by Archibald et al. (2020), these negative biases can be attributed to insufficient secondary production of CO from non-methane volatile organic compound oxidation and strong loss through hydroxyl radicals (OH) in the Northern Hemisphere. Archibald et al. (2020) also showed that, for UKESM1, there is a positive bias associated with regions where CO emissions are dominated by agricultural (eastern Central Asia region) and forest fires in central Africa (NHAF and SHAF) and the northwestern part of South America (NHSA). In UKESM1+INFERNO, the bias found in the CO column is similar to the ones found in UKESM1. However, the overestimation of biomass burning emissions of CO on the southern edge of NHAF and SHAF, previously described in Sect. 3.1, is reflected in a higher column CO positive bias in UKESM1+INFERNO when compared with TES-AURA for these regions. For South America (NHSA and SHSA regions), UKESM1+INFERNO shows an improvement in the western and central parts of NHSA but an increase of the negative bias in the southern and eastern parts of SHSA.

Figure 10 Column volume mixing ratio of carbon monoxide (ppm) averaged between 700 and 300 hPa for the period 2007–2012 (a) monthly mean time series and (b) monthly mean climatology for the NHAF region, (c) monthly mean time series and (d) monthly mean climatology for the SHAF region, (e) monthly mean time series and (f) climatology for the NHSA region, and (g) monthly mean time series and (h) monthly mean climatology for the SHSA region – UKESM1+INFERNO (solid line), UKESM1 (dotted line) and TES-AURA (dashed line).

A comparison of the monthly mean time series and monthly mean climatology between TES-AURA, UKESM1 and UKESM1+INFERNO for the main fire regions – NHAF, SHAF, NHSA and SHSA – is shown in Fig. 10 and provides information on interannual variability and the climatology of the regions that dominate the differences in the CO atmospheric column between UKESM1 and UKESM1+INFERNO. In most of these regions, the UKESM1+INFERNO configuration tends to provide an improvement on the modelled volume mixing ratio of CO, for NHAF (Fig. 10a and b), reducing the negative bias and RMSE of UKESM1 from −17.27 to −3.33 and 19.03 to 11.78 ppb, respectively. However, both configurations show a similar correlation of 96.28 in UKESM1 % and 95.92 % in UKESM1+INFERNO. These improvements are mostly associated with the decreases of the CO atmospheric column between May and August for UKESM1+INFERNO and the improvement of the variability over this region. As can be seen in Fig. 10c, SHAF is the region with the largest changes when comparing UKESM1 and UKESM1+INFERNO. When compared to TES-AURA, there is a significant worsening of the bias caused by the higher emissions over this region – from −8.46 to 20.70 ppb for UKESM1 and UKESM1+INFERNO, respectively. In addition, the RMSE is also increased from 19.03 ppb in UKESM1 to 28.88 ppb in UKESM1+INFERNO. Despite this, both configurations show similar high correlations with TES-AURA (94.89 % and 94.89 % for UKESM1 and UKESM1+INFERNO, respectively). These results show that, despite the large bias caused by the vegetation errors for this region, the UKESM1+INFERNO configuration captures the observed variability of CO, representing the two peaks in CO that occur in April and August. This can be seen both in monthly mean time series, as well as in the climatology for SHAF (Fig. 10d).

With regards to South America, there is an improvement in the NHSA region (Fig. 10e and f), with UKESM1+INFERNO presenting lower values of RMSE and bias than UKESM1 when compared to TES-AURA – RMSE reduction from 11.81 to 9.02 ppb and bias from −7.53 to −0.28 ppb – and the correlation presenting similar values 96.13 % for UKESM1 and 96.68 % for UKESM1+INFERNO. Even so, as depicted in Fig. 10f, for this region, UKESM1+INFERNO produces a more pronounced bimodal seasonality of CO atmospheric column peaking in late April and October, consistent with observations, but underestimates the magnitude especially for October, failing to represent years when concentrations are kept high during the summer months. It is worth noting that these results may be affected by the advection of CO from the SHAF region. In the SHSA region (Fig. 10g and h), UKESM1+INFERNO does not perform as well as UKESM1. Despite a good correlation with TES-AURA (91.94 %), the peak in CO volume mixing ratio tends to extend for longer periods of time throughout the year (from May to November) which results in a larger positive bias (8.55 ppb in UKESM1+INFERNO and −2.22 ppb in UKESM1). The results of the statistical comparison of the monthly mean time series, between the UKESM1 and UKESM1+INFERNO configurations with TES-AURA, are summarized in Table 3.

Table 3 Statistical comparison – root mean squared error (RMSE), bias and Pearson correlation coefficient (CORR) – of the monthly average time series between UKESM1 and UKESM1+INFERNO with TES-AURA for the period 2007–2012. Values in bold show that the correlation between a given experiment and TES-AURA is statistically significant with a 95 % confidence level.

3.4 Aerosols

Aerosols have a large effect on the Earth's radiative budget and climate; they can scatter and absorb radiation, as well as change cloud proprieties leading to changes in cloud cover and precipitation (Ward et al., 2012). As shown by Liousse et al. (1996), biomass burning is the primary source of natural carbonaceous aerosols in the ES (OC and BC), making them a major influence in controlling the ES variability in regions where fire activity is dominant. To evaluate the model's ability to reproduce the observed distribution and variability of aerosols, we compare the model to aerosol products from three instruments: MODIS and AERONET for AOD and CALIPSO for the vertical profile of the aerosol extinction coefficient. The focus is on the global distribution, seasonality and interannual variability of the major fire regions – NHAF, SHAF, NHSA and SHSA.

Mulcahy et al. (2019) provides a comprehensive model evaluation of aerosols in UKESM1 using prescribed fire emissions. In their study, the authors showed that the model performs well when compared to observations, capturing the global spatial distributions of AOD and cloud droplet number concentrations. The authors also report regional biases, including an overestimation of droplet number concentrations in the marine stratocumulus cloud regimes and an underestimation of aerosol optical depth in dust-dominated regions (Fig. 11c).

Figure 11 Aerosol optical depth multi-year average annual mean at 550 nm (2003–2012) for (a) MODIS, (b) UKESM1+INFERNO, (c) difference between UKESM1 and MODIS, and (d) difference between UKESM1+INFERNO and MODIS. Differences are only shown for statistically significant points with a 95 % confidence level.

Figure 12 Aerosol optical depth at 550 nm for the period 2003–2012 (a) monthly mean time series and (b) climatology of monthly means for the NHAF region, (c) monthly mean time series and (d) climatology of monthly means for the SHAF region, (e) monthly mean time series and (f) climatology of monthly means for the NHSA region, and (g) monthly mean time series and (h) climatology of monthly means for the SHSA region – UKESM1+INFERNO (solid line), UKESM1 (dotted line) and MODIS (dashed line).

This section focuses on the evaluation of UKESM1+INFERNO. When compared to MODIS, the annual mean of aerosol optical depth simulated by UKESM1+INFERNO (Fig. 11) has a realistic global spatial pattern. The global RMSE is 0.07, the bias is −0.07, and the global pattern correlation between the two model simulations is 79.7 %. Focusing on the main fire regions, the UKESM1+INFERNO tends to overestimate the AOD over Africa (both NHAF and SHAF), as well as SHSA and underestimate it in NHSA. Looking at the monthly mean time series and climatology of monthly means for AOD at 550 nm (Fig. 12), overall for the NHAF region (Fig. 12a and b), the UKESM1 and UKESM1+INFERNO results are dominated by a negative bias when compared to MODIS. This is found to be associated with the dust aerosol bias, as described in Mulcahy et al. (2019). Nonetheless, the southern edge of this region is dominated by fire emissions and both model configurations tend to overestimate the AOD, with an increase of the bias for UKESM1+INFERNO. Comparing the time series for this region (Fig. 12a), the higher biomass burning emissions of OC and BC in UKESM1+INFERNO compensates for the lack of dust emissions which results in a smaller bias and larger RMSE for UKESM1+INFERNO when compared to UKESM1 (bias of −0.06 from −0.10 and RMSE of 0.14 from 0.12). The pattern correlation of mean monthly values in UKESM1+INFERNO is also improved, when we compare to UKESM1 (from 21.02 % to 37.2 %). On the other hand, the AOD in the SHAF region (Fig. 12c and d) is dominated by biomass burning emissions and, as seen for CO, UKESM1+INFERNO tends to overestimate the biomass burning emissions due to the overestimation of tree fraction in this region. This leads to more OC and BC emissions. When comparing the time series between the two model configurations and MODIS (Fig. 12c), there is a change of signal, as well as an increase in the bias (from −0.003 in UKESM1 to 0.056 in UKESM1+INFERNO) and RMSE (from 0.05 for UKESM1 to 0.11 UKESM1+INFERNO). In terms of spatial correlation (Fig. 11), both models are highly correlated with MODIS observations, with UKESM1 having a higher correlation coefficient (84.3 %) than UKESM1 (71.9 %).

For both regions in South America, both model configurations present a very similar behaviour. In NHSA (Fig. 12e and f), both models reproduce well the AOD seasonal cycle (Fig. 12f), despite a constant bias that is present throughout the whole period. The interannual variability in both configurations is similar, with UKESM1+INFERNO performing slightly better compared to UKESM1, with a RMSE of 0.06 and a bias of −0.05 for UKESM1 and an RMSE of 0.05 and a bias of −0.03 for UKESM1+INFERNO. However, UKESM1+INFERNO does not reproduce some specific observed fire events that happen outside the normal fire season which are prescribed in UKESM1. Nonetheless, UKESM1+INFERNO has a better correlation than UKESM1 when compared to observations (60.5 % cf. 66.0 %). On a similar note, in SHSA (Fig. 12g and h), UKESM1 shows a RMSE of 0.05, a bias of −0.02 and a correlation of 88.3 %, while UKESM1+INFERNO shows a RMSE of 0.06, a bias of −0.006 and a correlation of 80.9 %. UKESM1+INFERNO does not reproduce the large AOD peaks that occur during the period 2004 to 2007 and in 2010 (Fig. 12g) – which are also present for CO (Fig. 10g). These peaks in AOD are associated with the Amazonian fire events that generally occur in drought years, which are often related to El Niño events. The ability of the model to reproduce these specific events depends on the ability to represent circulation regimes that lead to the drought but, more importantly, the fire ignitions associated human activities other than deforestation including secondary vegetation slash-and-burn and cyclical fire-based pasture cleaning (which are boosted in drought years) which are not represented in INFERNO (Aragão et al., 2018; Marengo et al., 2011).

The AERONET Sun photometers provide a ground-based direct measurement of the attenuation of sunlight due to aerosol. This means they are not affected by the same uncertainties as satellite retrievals, associated with the different satellite retrieval algorithms (e.g. assumptions related to underlying surface properties as in MODIS). Despite some AERONET stations providing the most extensive records of AOD observations, the location of these stations is sparse in many of the key regions where aerosols are dominated by biomass burning emissions. For this reason, a compromise between a long record and including stations within the regions of interest for this analysis was found, and all AERONET sites that have at least a 5-year continuous record that overlap with model data were included.

Figure 13 Aerosol optical depth mean seasonal average at 440 nm (1980–2014) for the climatological seasons of DJF, MAM, JJA and SON (represented in each row) for UKESM1 (a, c, e, g) and UKESM1+INFERNO (b, d, f, h). The ground-based aerosol optical depth retrievals at various AERONET sites that have at least a 5-year overlapping period with model data are overlaid in circles using the same colour scale.

When comparing modelled AOD results with AERONET (Fig. 13), it is possible to see that these agree broadly with the analysis previously done for MODIS AOD. However, contrary to MODIS where UKESM1+INFERNO shows a negative bias, and possibly related to MODIS retrieval bias over bright surfaces, stations located in NHAF suggest that the model performs well during the boreal winter (DJF) (Fig. 13b) and underestimates the AOD for spring (MAM) (Fig. 13d). Moreover, for this region, there is a higher AOD during MAM for UKESM1+INFERNO (Fig. 13d) than for UKESM1 (Fig. 13c), bringing this model configuration closer to the AOD values found for the nearby AERONET station. On a similar note, for the SHAF region, there is also a better agreement between UKESM1+INFERNO and the AERONET AOD during the SON season (Fig. 13h), with UKESM1 (Fig. 13g) showing a negative bias relative to the stations found for this region. As seen before for South America (both NHSA and SHSA), both model configurations exhibit similar performance and show similar biases when comparing to MODIS – an overall underestimation of the AOD over these regions during the biomass burning season (SON).

Further to the differences introduced by the INFERNO interactive biomass burning emissions in UKESM1+INFERNO described previously, a different approach to the one adopted in UKESM1 has been taken regarding the vertical distribution of these emissions. For consistency between the approach regarding gas-phase biomass burning emissions, in UKESM1+INFERNO, these are injected into the model's lowest layer. However, while in UKESM1 prescribed aerosol biomass burning emissions are distributed uniformly over the 20 first model levels (∼3 km), in UKESM1+INFERNO these emissions are distributed following an exponential increasing with height function from the first model level to the 20th model level. In order to evaluate the vertical profile of modelled aerosol, the CALIPSO lidar level 3 aerosol extinction coefficient profiles product at 532 nm were used. Following the analysis of previous sections, we focused this analysis on the African (NHAF and SHAF) and South American (NHSA and SHSA) regions.

Figure 14 Aerosol extinction coefficient mean (km⁻¹) annual average vertical cross section at 532 nm (2006–2014) covering the African continent (extending both NHAF and SHAF regions) for (a) CALIPSO, (b) UKESM1+INFERNO, (c) difference between UKESM1 and CALIPSO, and (d) difference between UKESM1+INFERNO and UKESM1. Stapling is shown for points where the differences are statistically significant with a 95 % confidence level.

Over Africa – Fig. 14a – CALIPSO shows larger extinction coefficient in the first kilometre of the atmosphere, extending vertically up to 5 km over the tropical regions (5–15^∘ N in latitude), associated with the vertical transport of aerosols. There are two main areas of large extinction coefficients: one associated with dust aerosol emissions from the Sahel region – down to ∼10^∘ N as well as the southern edge biomass burning emission of NHAF extending further down to 0^∘ N – and a second area from −5^∘ N extending south to −25^∘ N associated with the SHAF fire region. When comparing the model to CALIPSO (Fig. 14c and d), it can be seen that both of these model configurations show identical differences. On the lower levels (for altitudes below 1 km), there are significant large negative biases (>75 % compared to CALIPSO), except for biomass burning area around 0^∘ N, where there is a small overestimation (<10 %). These lower-level biases also affect higher levels as aerosol are transported vertically. However, for altitudes above 1 km over the regions of the large negative bias, there is a reduction to 25 %–50 % when compared to CALIPSO. In addition, for latitudes below −25^∘ N, the models result show a positive bias (<0.05 km⁻¹), extending from 1 to 5 km in the vertical.

When comparing the extinction coefficient vertical profile of UKESM1+INFERNO to UKESM1 – Fig. 14d – it can be seen that the differences are 1 order of magnitude smaller than those found when comparing UKESM1+INFERNO to CALIPSO. UKESM1+INFERNO shows larger values of extinction coefficient at higher levels (above 1.5 km) and smaller at the lower levels (below 1.5 km) when compared to UKESM1 (Fig. 14c); this could be both an impact of the different treatment of emissions in the vertical, as well as due to transport associated with different locations for the biomass burning emissions. In addition, and as discussed before, due to the further extension towards south of the SHAF fire region, there is a positive bias throughout the vertical levels for latitudes southward of −20^∘ N.

Figure 15 Aerosol extinction coefficient mean annual average vertical cross section at 532 nm (2006–2014) covering the South American continent (extending both NHSA and SHSA regions) for (a) CALIPSO, (b) UKESM1+INFERNO, (c) difference between UKESM1 and CALIPSO, and (d) difference between UKESM1+INFERNO and UKESM1. Stapling is shown for points where the differences are statistically significant with a 95 % confidence level.

Similar results can be found when comparing UKESM1+INFERNO to CALIPSO for South America (Fig. 15). South American aerosol emissions are dominated by biomass burning, especially around tropical regions (−15 to 15^∘ N). There is a large underestimation of the extinction coefficient for altitudes below 1 km (bias between 25 % and 50 % of the observed value). However, contrary to what could be seen for Africa, there is a positive bias for altitudes above 1 km (<0.02 km⁻¹) in this region. When comparing UKESM1+INFERNO and UKESM1 – Fig. 15d – the difference is dominated by a large negative bias (∼25 % of the observed value in CALIPSO, as can be seen in Fig. 15a) between −25 and −5^∘ N extending from the surface to an altitude of 4 km. This difference is consistent with the results found for AOD over this region, with UKESM1 (Fig. 15c) presenting a higher AOD when compared to UKESM1+INFERNO. This result, together with the differences between UKESM1+INFERNO and UKESM1 found for Africa, suggest that the influence of the underlying model bias in aerosol vertical and spatial distributions is dominant when compared to the influence of the different treatments of the vertical distribution of biomass burning aerosols.

4 Discussion and conclusions

The goal of this work was the development and evaluation of the implementation of a coupled fire–climate–composition ES model. This was built on top of the work developed by Mangeon et al. (2016), coupling the INFERNO fire model to the atmosphere-only configuration of version 1 of the UK's Earth System Model (UKESM1). The fire–atmosphere interactions through atmospheric chemistry and aerosols allow for fire emissions to provide feedback on radiation and clouds, changing weather, which can consequently provide feedback on the atmospheric drivers of fire. This is the basis for the development of a framework that would allow the impacts of fire variability on atmospheric composition–climate interactions in the past, present and future in an ES model context to be quantified. It also provides the possibility of adding further coupling to a dynamic land surface vegetation model in order to capture the fire–vegetation–atmospheric composition–climate interactions.

During the development of this work, it was identified that INFERNO human fire ignitions and fire suppression functions excluded the representation of socio-economic factors (aside from population) that can affect anthropogenic behaviour regarding fire ignitions. To address this, we include an HDI term aimed at representing socio-economic factors impacting fire ignition and suppression. The HDI is calculated based on three indicators designed to capture the income, health and education dimensions of human development. Therefore, we assume this leads to a representation where if there is more effort in improving human development, there is also investment on higher fire suppression by the population. Furthermore, the biomass burning emissions factors were updated according to the work of Andreae (2019) to reflect the state-of-the-art knowledge available since the development of INFERNO and emissions for C₂H₆, C₃H₈, HCHO, MeCHO, Me2CO and DMS were also added to make the fire–composition coupling consistent between UKESM1 and UKESM1+INFERNO. One of the limitations of the current biomass burning emission coupling framework is that the vertical distribution of emissions is not dependent on the radiative power of fire events. Instead, all fire events apply the same treatment to distribute these emissions vertically in the atmospheric column – surface emissions for gas-phase emissions and a linear increasing profile for aerosol emissions. Nevertheless, the results presented here suggest that the background bias in aerosols outweighs the differences introduced by using a different emission profile (Figs. 14 and 15). Moreover, it was shown that INFERNO is highly sensitive to the underlying land surface types of the land surface provided to the model. The prescribed land surface types used in this configuration of UKESM1+INFERNO overestimates the tree fraction in savanna biomes, such as the southern Africa (SHAF) and the southern edge of the Amazon rain forest (SHSA), which decreased the potential average burnt area and increased the fuel available for combustion and emissions, causing notable changes to the fire seasonality. This highlights the importance of correctly simulating land cover and it suggests that including the coupling between fire and vegetation in this framework could significantly improve model results as it would allow fires to shape regional land cover.

The model has a realistic spatial distribution of the average burnt area at a continental scale. However, there is a large underestimation of the annual average burnt area of approximately 250 Mha due to the underestimation of fires in NHAF and SHAF caused by the overestimates in the tree fraction for these regions. This underestimation of burnt areas impacts the biomass burning emissions. Despite the overall global pattern for biomass burning emissions being well reproduced by the model, there is a large overestimation of the biomass burning emissions for all the emitted species on the southern edge of NHAF and SHAF, as well as the eastern SHSA side (difference of > 300 %). Furthermore, in the SHAF region, the emissions extend further south into the midlatitudes. At a global scale, the effects of the underestimation of burnt areas is compensated by the increased availability of carbon provided by the overestimation of the tree fraction, resulting in good agreement with the observations.

Comparing UKESM1+INFERNO to TES-AURA satellite CO product showed that, despite the overestimation of CO over NHAF and SHAF due to the bias in the emissions, UKESM1+INFERNO has a similar performance when compared to UKESM1. In fact, including the interactive biomass burning emissions improves the interannual CO atmospheric column variability and consequently its seasonality over the main biomass burning regions – Africa and South America. Similarly, for aerosols, the AOD results broadly agree with MODIS and AERONET observations. Most of the biases found for aerosols are also present in UKESM1 and are not associated with the interactive biomass burning emissions. When comparing UKESM1+INFERNO to the observed datasets, and taking into account the background bias found in UKESM1, it can be seen that there is an overestimation of AOD over the biomass burning regions in Africa and an improvement of the variability and seasonality of AOD in South America, with UKESM1+INFERNO presenting better correlation and lower bias than UKESM1 for SHSA. Nonetheless, UKESM1+INFERNO does not capture the spikes in AOD, or CO observed over SHSA during the period 2004 to 2007 and 2010 associated with the Amazonian fire events that generally occur in drought years often related to El Niño events. Capturing these specific events depends on the ability of the model to represent circulation regimes that lead to the drought during the simulated period, as well as the vegetation flammability and fire ignitions. When analysing the vertical profile of aerosol extinction coefficient comparing to CALIPSO data, it was possible to see that there is an underestimate of aerosols in the first kilometre of the atmosphere which is present both in UKESM1 and UKESM1+INFERNO. This suggests that this is an underlying bias of the UKESM1 configuration – not related to the coupling of biomass burning emissions. These results are consistent with the results found for AOD over these regions, where UKESM1 is found to have a higher AOD when compared to UKESM1+INFERNO. This result, together with the differences between UKESM1+INFERNO and UKESM1 found for Africa suggest that the underlying model biases regarding aerosol spatial distributions are dominant when compared to the different treatment of the vertical distribution of the biomass burning aerosols.

Despite the present limitations of the implementation of coupled fire–composition–climate processes in this framework, it is evident that the UKESM1+INFERNO demonstrates a similar performance in reproducing the distribution of aerosols and CO atmospheric column. This shows that UKESM1+INFERNO provides a useful coupling framework that allows an internally consistent representation of complex fire–atmospheric composition–climate complex interactions and feedbacks in the Earth system. With further work to improve the existent model errors and bias, such as the coupling the fire feedbacks to vegetation and representing peatland regions, UKESM1+INFERNO can provide a useful framework to quantify the impacts of fire variability on atmospheric composition–climate interactions in past, present and future climates.