Studying the impact of biomass burning aerosol radiative and climate effects on the Amazon rainforest productivity with an Earth system model

Abstract

Diffuse light conditions can increase the efficiency of photosynthesis and carbon uptake by vegetation canopies. The diffuse fraction of photosynthetically active radiation (PAR) can be affected by either a change in the atmospheric aerosol burden and/or a change in cloudiness. During the dry season, a hotspot of biomass burning on the edges of the Amazon rainforest emits a complex mixture of aerosols and their precursors and climate-active trace gases (e.g. ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{CH}}_{\mathrm{4}}$, ${\mathrm{NO}}_x$). This creates potential for significant interactions between chemistry, aerosol, cloud, radiation and the biosphere across the Amazon region. The combined effects of biomass burning on the terrestrial carbon cycle for the present day are potentially large, yet poorly quantified. Here, we quantify such effects using the Met Office Hadley Centre Earth system model HadGEM2-ES, which provides a fully coupled framework with interactive aerosol, radiative transfer, dynamic vegetation, atmospheric chemistry and biogenic volatile organic compound emission components. Results show that for present day, defined as year 2000 climate, the overall net impact of biomass burning aerosols is to increase net primary productivity (NPP) by $+\mathrm{80}$ to $+\mathrm{105}$ TgC yr${}^{-\mathrm{1}}$, or 1.9 % to 2.7 %, over the central Amazon Basin on annual mean. For the first time we show that this enhancement is the net result of multiple competing effects: an increase in diffuse light which stimulates photosynthetic activity in the shaded part of the canopy ($+$65 to $+$110 TgC yr${}^{-\mathrm{1}}$), a reduction in the total amount of radiation ($-\mathrm{52}$ to $-\mathrm{105}$ TgC yr${}^{-\mathrm{1}}$) which reduces photosynthesis and feedback from climate adjustments in response to the aerosol forcing which increases the efficiency of biochemical processes ($+$67 to $+$100 TgC yr${}^{-\mathrm{1}}$). These results illustrate that despite a modest direct aerosol effect (the sum of the first two counteracting mechanisms), the overall net impact of biomass burning aerosols on vegetation is sizeable when indirect climate feedbacks are considered. We demonstrate that capturing the net impact of aerosols on vegetation should be assessed considering the system-wide behaviour.

Introduction

The Amazon rainforest is the largest expanse of tropical forest on Earth. It provides invaluable ecological services and plays a major role in the Earth system and climate (Malhi et al., 2008). The Amazon rainforest is a net sink of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$, although drought frequency and intensity, which are expected to increase in the future, could have severe consequences for future forest resilience and potentially shift the Amazon rainforest from a sink to a net source of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ (Cox et al., 2000, 2004; Phillips et al., 2009; Doughty et al., 2015; Duffy et al., 2015; Sakschewski et al., 2016; Zemp et al., 2017). This possibility motivated intense research to develop a better understanding of the rainforest response to environmental stresses via integrated explicit representations of the carbon cycle in Earth system models (ESMs) (Cox et al., 2000). Response to many of these environmental stresses is now well documented and represented in ESMs, including the effects of surface temperature, atmospheric composition, water availability, or the amount and quality of accessible light (direct versus diffuse) for plant photosynthesis (e.g. Nemani et al., 2003; Sitch et al., 2007; Cox et al., 2008; Mercado et al., 2009a; Beer et al., 2010; Ciais et al., 2013; Pacifico et al., 2015; Unger et al., 2017).

In parallel to the above-mentioned environmental stresses, forest fires are also an intrinsic component of some forest lifecycles, providing an additional mechanism for depleting land carbon reservoirs. Intense biomass burning events present a notorious pressure on tropical regions and typically occur during the dry season – i.e. between around August and September in the Amazon region (Artaxo et al., 2013; Brito et al., 2014). Fires in general occur naturally; however, a significant fraction results from the anthropogenic pressure that continually erodes the fragmented forest edges (Cochrane, 2003). Despite a decreasing trend in the rate of deforestation over the last decade as a result of stricter environmental policies (Kalamandeen et al., 2018), it is estimated that 293  Tg of carbon per year (TgC yr${}^{-\mathrm{1}}$) is directly released back into the atmosphere from fires in the Amazon (van der Werf et al., 2006). Fires can also have an indirect impact on the rainforest carbon budget that is harder to quantify; for instance, fires alter surface properties (e.g. albedo) in the burnt area, which can modify surface fluxes and the water cycle (e.g. Zemp et al., 2017). Additionally, fires emit a complex mixture of gases (${\mathrm{CO}}_{\mathrm{2}}$, $\mathrm{CO}$, ${\mathrm{CH}}_{\mathrm{4}}$, ${\mathrm{NO}}_x$ and volatile organic compounds – VOCs), aerosols and aerosol precursors which can affect remote regions of the rainforest after being dispersed by the wind. Pacifico et al. (2015) illustrated such a mechanism by analysing the potentially harmful effect of near-surface ozone (${\mathrm{O}}_{\mathrm{3}}$) associated with biomass burning and estimated that the rainforest gross primary productivity (GPP) was reduced by up to approximately $-\mathrm{230}$ TgC yr${}^{-\mathrm{1}}$, a number of similar magnitude to the direct carbon loss from fires.

Assessing the overall impact of Amazonian forest fires on ecosystems is challenging as it encompasses a combination of direct losses and indirect impacts from the fire by-products which can depend on intricate interactions among several Earth system components, including the biosphere, atmospheric composition, radiation and energy budget, clouds, and the water cycle (Bonan, 2008). Here, we aim to specifically elucidate the impact of biomass burning aerosols (BBAs) that are associated with forest fires and quantify their potential effect on the Amazon forest productivity.

Significant amounts of BBAs are emitted in South America, which strongly modify the radiative budget by scattering and absorbing solar radiation. This reduces the level of photosynthetically active radiation (PAR), traditionally defined as the radiation between wavelengths of 300 and 700 nm, reaching the surface and used by plants to photosynthesise (i.e. to assimilate carbon from the atmosphere). Contrary to intuition, an increase in the diffuse light fraction can be beneficial to plants as the shaded, non-light-saturated leaves, typically found in the understory or lower canopy layers, receive more radiation under diffuse light conditions than they would normally experience under direct light conditions owing to the shading by leaves fully exposed to sunlight. As a result, this trade-off between experiencing less PAR overall and receiving more evenly distributed light across the canopy favours higher rates of canopy photosynthesis. The first comprehensive estimation of this diffuse PAR fertilisation effect (DFE) at the global scale was documented by Mercado et al. (2009a), who used a combination of offline aerosol distributions, radiative transfer and a land surface model to estimate that DFE may have increased the global land carbon uptake by up to 25 % during the global dimming period (1950–1980; Stanhill and Cohen, 2001). More recently, Rap et al. (2015) used a similar framework of offline models to assess the role of BBA over the Amazon region. They showed that BBA increases the annual mean diffuse light and net primary productivity (NPP) by 3.4 %–6.8 % and 1.4 %–2.8 %, respectively. Strada and Unger (2016) took a step further using a coupled modelling framework to estimate biomass burning aerosol impacts on Amazon forest GPP, obtaining an increase of 2 %–5 % on annual means. Recently, Moreira et al. (2017) also applied a coupled framework using a regional model (BRAMS) to conclude that BBA could increase the GPP of the Amazon forest by up to 27 % during the peak of the biomass burning season. The study of Moreira et al. (2017) assumed high BBA emissions and did not account for the effect of cloudiness on the diffuse fraction of radiation, so it provides an upper estimate of the potential impact of the effects of the attenuation of total solar radiation and the enhancement of the diffuse solar radiation flux inside the vegetation canopy.

Despite a growing body of evidence supporting the DFE mechanism, both from observational and modelling perspectives (e.g. Cohan et al., 2002; Gu et al., 2003; Robock et al., 2005; Yamasoe et al., 2006; Mercado et al., 2009a; Kanniah et al., 2012; Cirino et al., 2014; Cheng et al., 2015), a full quantification of the BBA impact on ecosystems remains poor because aerosol–radiation interactions (ARI), and to some extent aerosol–cloud interactions (ACI), not only create the conditions for a DFE but also modify the climate locally. For example, a regional haze of aerosols can perturb regional hydroclimates (Nigam and Bollasina, 2010), force clouds to adjust to aerosol semi-direct and indirect effects which modify the way clouds interact with radiation (Hansen et al., 1997; Haywood and Boucher, 2000; Koren et al., 2004), or create a positive cooling effect on productivity by reducing surface heat stress in hot environments, allowing for a more efficient uptake of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ through leaf stomata (Robock et al., 2005; Xia et al., 2016; Strada and Unger, 2016). Neglecting such essential coupling pathways may overemphasise the relative contribution of the DFE due to loss of internal consistency that does not allow variability within non-linear relationships. Only a limited number of studies have considered the DFE within a fully coupled Earth system framework (e.g. Strada and Unger, 2016; Unger et al., 2017; Yue et al., 2017, using the NASA GISS ModelE2–YIBs) to investigate the role of aerosols and haze on vegetation. Although these studies have investigated the role of diffuse radiation on GPP and isoprene emissions (Strada and Unger, 2016; Unger et al., 2017), understanding of the indirect impact of climate effects from aerosols on vegetation productivity remains very uncertain. This was addressed over China by Yue et al. (2017), who demonstrated that aerosol-induced hydroclimatic feedbacks can promote ecosystem NPP. In the present study, we apply an ESM modelling framework to quantify the impact of present-day BBA via the quantification of individual and net effects of changes in diffuse radiation, direct radiation and climate upon the vegetation productivity in the Amazon rainforest specifically. For this endeavour, we have implemented an updated representation of plant photosynthesis and carbon uptake that is sensitive to diffuse light radiation in the UK Met Office Hadley Centre HadGEM2-ES Earth system model (Mercado et al., 2007, 2009a). In addition, a framework that disentangles the vegetation response has been developed to provide a deeper understanding of the contributions of different plant environmental variables affected by aerosols. The role of ${\mathrm{O}}_{\mathrm{3}}$ precursor emissions and in situ formation of ${\mathrm{O}}_{\mathrm{3}}$ associated with biomass burning (Pacifico et al., 2015) is not considered here.

The methodology and the experimental set-up are described in Sect. 2. Results are discussed in Sect. 3, including first a model evaluation in Sect. 3.1, then the net effect of BBA in Sect. 3.2, and individual contributions from the diffuse light fraction, the reduction in total PAR and the climate feedbacks associated with the BBA perturbation in Sect. 3.3. These findings are contextualised in Sect. 3.4 by analysing the results from four additional sensitivity experiments designed to elucidate the role of aerosol optical properties, aerosol–cloud interactions, the atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration and vertical distribution of nitrogen through the canopy. Concluding remarks and a summary of this study's main results are provided in Sects. 4 and 5, respectively.

Method

We evaluate the effects of biomass burning aerosol–radiation interactions upon the Amazon rainforest primary productivity for present-day conditions using the Met Office Hadley Centre Global Environment Model HadGEM2-ES (The HadGEM2 Development Team, 2011), which provides a fully coupled framework. The model is briefly described in Sect. 2.1.

We present the results of a sensitivity experiment (Sect. 3) which consists of varying the biomass burning aerosol emissions only over South America. “Real world” fires also emit greenhouse gases (e.g. ${\mathrm{CO}}_{\mathrm{2}}$, $\mathrm{CO}$, ${\mathrm{CH}}_{\mathrm{4}}$) and ozone precursors (${\mathrm{NO}}_x$, VOCs) which can potentially affect the biosphere. Ozone is particularly critical as it is a pollutant which harms plants and reduces their productivity, and thus their ability to draw ${\mathrm{CO}}_{\mathrm{2}}$ from the atmosphere (Sitch et al., 2007). Whereas the damaging effect of ozone is not accounted for in this study, we will briefly discuss the potential fertilisation effect from the increased ${\mathrm{CO}}_{\mathrm{2}}$ background that can result from biomass burning in Sect. 4. The ozone damage effect was documented by Pacifico et al. (2015) using a similar modelling framework as in the present study, and we refer readers to that study for further details.

Atmospheric particles such as aerosols and cloud droplets scatter radiation, which increases the fraction of radiation that is diffuse. Diffuse conditions result in higher light use efficiency of plant canopies, which can enhance carbon uptake (Roderick et al., 2001; Gu et al., 2002). An increase in diffuse radiation is concomitant with a decrease in the overall amount of radiation (Fig. S1 in the Supplement). These two opposing effects will be referred to in the rest of the paper as “change in diffuse fraction” and “reduction in total PAR”, respectively, and will be quantified separately in Sect. 3.3. Finally, BBA effects impact the coupled system, which controls the rate of biochemical processes of vegetated land surfaces itself. We will simply refer to these adjustments to the BBA effects as “climate feedback” in the remainder of the paper. The sum of climate feedback, change in diffuse fraction and reduction in total PAR is referred as the “net impact” of BBA on plant productivity. The framework we developed to disentangle these three terms is described in Sect. 2.4.

Model description

HadGEM2-ES is an Earth system model built around the HadGEM2 atmosphere–ocean general circulation model and includes a number of Earth system components such as

• the ocean biosphere Diat-HadOCC (Diatom-Hadley Centre Ocean Carbon Cycle) model, developed from the HadOCC model of Palmer and Totterdell (2001);

• the sea ice component (The HadGEM2 Development Team, 2011);

• the Top-down Representation of Interactive Foliage and Flora Including Dynamics (TRIFFID) dynamic global vegetation model (Cox, 2001), and the land surface and carbon cycle model MOSES2 (Met Office Surface Exchange Scheme), collectively known as JULES (Cox et al., 1998, 1999; Essery et al., 2003);

• the interactive Biogenic Volatile Organic Compounds (iBVOC) emission model (Pacifico et al., 2012);

• the UKCA tropospheric chemistry scheme (O'Connor et al., 2014).

The atmospheric model resolution is N96 (1.875${}^{\circ }$ by 1.25${}^{\circ }$) with 38 vertical levels with the model top at $\sim \mathrm{39}$ km. Our modelling framework is similar to the configuration used by Pacifico et al. (2015), who provided a detailed analysis of the successful model performance against observations.

For clarity, we provide some additional details on the treatment of aerosols and their coupling with radiation and clouds as well as on the updated representation of the canopy interaction with radiation. The radiative transfer code in the atmospheric part of HadGEM2-ES is SOCRATES (Edwards and Slingo, 1996), which parameterises radiative fluxes using a “two-stream” approximation (Meador and Weaver, 1980). The radiative transfer is solved for six wavebands in the shortwave and nine in the longwave. This scheme accounts for the interaction of radiation with aerosol particles by defining three single scattering properties on a layer: optical depth, single scattering albedo (the ratio of scattering efficiency to total extinction) and an asymmetry parameter. Together, these properties determine the overall transmission and reflection coefficients of each atmospheric layer. At the interface between the lowest atmospheric level and the land surface, the total and the direct radiances for the shortwave band 320–690 nm, which approximates the PAR, calculated by the SOCRATES radiation scheme are transferred to the land surface routines to calculate plant photosynthesis.

In the JULES land surface model, the total and direct irradiance components of PAR calculated by the atmospheric model provide the boundary conditions at the top of the canopy. The diffuse PAR fraction is calculated as the difference between total and direct radiation, divided by the total radiation. The canopy is discretized into 10 vertical layers, and the radiative transfer in the canopy is also parameterised with a two-stream approximation but uses more detailed assumptions to represent light interception by foliage (Sellers, 1985). The photosynthesis model is based upon the observed processes of gas and energy exchange at the leaf scale, which are then scaled up to represent the canopy. It takes into account variations in direct and diffuse radiation on sunlit and shaded canopy photosynthesis at each canopy layer. In this way, photosynthesis of sunlit and shaded leaves is calculated separately under the assumption that shaded leaves receive only diffuse light and sunlit leaves receive both diffuse and direct radiation (Dai et al., 2004; Clark et al., 2011). Leaf-level photosynthesis is calculated using the biochemistry of C3 and C4 photosynthesis from Collatz et al. (1991, 1992).

This canopy radiation scheme was first developed to quantify the impact of anthropogenic aerosol emissions on the global carbon cycle (Mercado et al., 2007, 2009a) and was consequently implemented in JULES (Clark et al., 2011). It is a novel addition to HadGEM2-ES as it was not available during the HadGEM2-ES contribution to CMIP5. HadGEM2-ES with the previous canopy radiation scheme had a tendency to overestimate GPP (Shao et al., 2013), which has to be balanced by high plant respiration (RESP) to get satisfactory estimates of global NPP (i.e. NPP $=$ GPP-RESP). The new representation of light interception that we have implemented is able to reproduce higher light use efficiency (LUE) under diffuse light conditions (Sect. 3.1 and Fig. S2 in the Supplement). However, the ratio of GPP to plant respiration in HadGEM2-ES with the new canopy radiation model remains too high when compared to observationally based estimates (e.g. Luyssaert et al., 2007). To correct this deficiency, we decreased the ratio of nitrogen allocated in the roots relative to the nitrogen in the leaves from 100 % to 50 % (Clark et al., 2011, Table 2 therein). Additionally, we reduced the leaf dark respiration coefficient that relates leaf dark respiration and $V_{\mathrm{cmax}}$ from 15 % to 10 % (Clark et al., 2011, Eq. 13 therein). These changes are based on a sensitivity analysis that we performed with the stand-alone version of JULES. We used the meteorological observations from the tropical French Guiana site (assumed to be fully covered by broadleaf trees) to drive JULES and investigate the sensitivity to parameters such as the leaf nitrogen content at canopy top (N${}_{\mathrm{L}\mathrm{0}}$), the dark respiration coefficient and the nitrogen allocation throughout the canopy via the value of the nitrogen profile extinction coefficient (Clark et al., 2011, Eq. 33 therein and Sect. 2.3.4 of the present study). Fast carbon fluxes (GPP, RESP and NPP) were calculated at a 3 h temporal resolution by varying one of these three parameters individually (Fig. S3a–c) and then averaged to annual mean values (Fig. S3d–f). The annual means were then used to construct contour surfaces for the fast carbon fluxes by varying combinations of the selected parameters (Fig. S4). This method enables us to ultimately pre-calibrate the fast carbon fluxes in the HadGEM2-ES model offline.

Aerosols are represented by the CLASSIC aerosol scheme (Bellouin et al., 2011) which is a one-moment mass prognostic scheme. This aerosol module contains numerical representation of up to eight tropospheric aerosol species. Here, ammonium sulfate, mineral dust, sea salt, fossil fuel black carbon (FFBC), fossil fuel organic carbon (FFOC), biomass burning aerosols and secondary organic (also called biogenic) aerosols are considered. Dust and sea salt are from diagnostic schemes based on the near-surface wind speed, while other emissions including biogenic aerosols are represented by a relatively simple climatology (Bellouin et al., 2011). Transported species experience boundary layer and convective mixing and are removed by dry and wet deposition. Wet deposition by large-scale precipitation is corrected for re-evaporation of precipitation: tracer mass is transferred from a dissolved mode to an accumulation mode in proportion to re-evaporated precipitation. For convective precipitation, accumulation mode aerosols are removed in proportion to the simulated convective mass flux. Emissions of biomass burning aerosols are the sum of the biomass burning emissions of black and organic carbon. Grass fire emissions are assumed to be located at the surface, while forest fire emissions are injected homogeneously across the boundary layer (0.8–2.9 km).

The direct radiative effect due to scattering and absorption of radiation by all eight aerosol species represented in the model is included. The semi-direct effect, whereby aerosol absorption tends to change cloud formation by warming the aerosol layer, is thereby included implicitly. Wavelength-dependent specific scattering and absorption coefficients are obtained using Mie calculations from prescribed size distributions and refractive indices. All aerosol species except mineral dust and fossil fuel black carbon are considered to be hydrophilic, act as cloud condensation nuclei, and contribute to both the first and second indirect effects on clouds, treating the aerosols as an external mixture. Jones et al. (2001) detailed the parameterization of the indirect effects used in HadGEM2-ES. The cloud droplet number concentration (CDNC) is calculated from the number concentration of the accumulation and dissolved modes of hygroscopic aerosols. For the first indirect effect, the radiation scheme uses the CDNC to obtain the cloud droplet effective radius. For the second indirect effects, the large-scale precipitation scheme uses the CDNC to compute the auto-conversion rate of cloud water to rainwater (Jones et al., 2001).

Experimental design: main experiment

The HadGEM2-ES model was initiated on 1 December 2000 from a previous historical simulation. We consider the year 2000 to be a good surrogate for present-day climate, which will enable us to assess the impact of present-day BBA emissions on vegetation. As historical simulations are transient climate simulations, we constrain the carbon cycle to present-day values as well (to be described in the next paragraph). The model is then integrated for a period of 40 years using periodic forcing for the year 2000 to construct an ensemble that captures the model internal variability. Results reported here are the multi-annual means over the final 30 years of the model integration. The domain of analysis is defined by the coordinates 0–15${}^{\circ }$ S, 70–53${}^{\circ }$ W and is primarily covered by broadleaf trees for this configuration of HadGEM2-ES (Fig. S5).

The HadGEM2-ES model is set-up in an Atmospheric Model Intercomparison Project (AMIP; Jones et al., 2011) type configuration using prescribed climatologies of monthly mean sea surface temperatures (SSTs) and sea ice cover (SIC), which enables us to analyse the rapid adjustments of land surface climate to aerosol radiation perturbations. The introduction of a new canopy radiation interaction model introduces a significant departure in the carbon cycle balance. To prevent the need of a complex spin-up exercise, we prescribe the vegetation cover and carbon reservoirs to present-day level. This is achieved by reducing the call frequency of the TRIFFID dynamic vegetation model to 30 years in order to maintain the vegetation in a steady state. A similar approach is discussed in Strada and Unger (2016). Overall, this enables us to focus our analysis on the fast carbon flux responses (i.e. NPP, GPP) and their sensitivity to the perturbation induced by the biomass burning aerosols.

Aerosols and their precursor emissions are the dataset used during CMIP5 (Lamarque et al., 2010). We use the decadal mean emissions centred around the year 2000 to represent present-day emission rates. Biogenic volatile organic compound (BVOC) emissions from vegetation (Pacifico et al., 2012) are sensitive to changes in plant productivity and hence sensitive to DFE. These emissions are calculated online but are not taken into account in the CLASSIC aerosol scheme. Instead, the climatology of BVOCs (also called secondary organics) from CMIP5 is used. The biomass burning emissions are based on the GFEDv2 inventory (van der Werf et al., 2006; Lamarque et al., 2010). Given the substantial inter-annual variability of biomass burning on a global and regional scale, a present-day climatology (i.e. average year) is calculated as the GFEDv2 1997–2006 average (Lamarque et al., 2010). These are the standard emission scenarios for the simulation labelled as BBAx1 for the main experiment. A total of five simulations are conducted in the main experiment where the standard biomass burning aerosols emissions are varied by $-\mathrm{100}$ %, $-\mathrm{50}$ %, 0 %, $+\mathrm{100}$ % and $+\mathrm{300}$ %, respectively (simulation BBAx0, BBAx0.5, BBAx1, BBAx2 and BBAx4, respectively). A multiplication factor is applied to the emission only for the BB sources over South America (40${}^{\circ }$ S, 85${}^{\circ }$ W; 15${}^{\circ }$ N, 30${}^{\circ }$ W). We define the control simulation as the simulation without BBA being emitted over South America (i.e. BBAx0). The changes in fast carbon fluxes are calculated as the departure from this reference simulation (e.g. $\mathrm{\Delta }$ NPP${}_{\mathrm{net}\mathrm{impact}}^{\mathrm{BBAx}\mathrm{1}}={\mathrm{NPP}}^{\mathrm{BBAx}\mathrm{1}}-{\mathrm{NPP}}^{\mathrm{BBAx}\mathrm{0}}$ and represents the net change in NPP due to standard emissions of BBA).

Sensitivity experiments

In parallel to the five simulations for the main experiment, we have conducted the following four additional sensitivity experiments to further appreciate the role of (i) aerosol optical properties, (ii) aerosol–cloud interactions, (iii) the canopy nitrogen profile and (iv) atmospheric carbon dioxide concentration. A listing of the simulations done for the main experiment and the sensitivity experiments is provided in Table 1.

List of model simulations done for the five experiments.

			Main	Sensitivity experiment:		Sensitivity experiment:		Sensitivity experiment:	Sensitivity experiment:
			experiment	aerosol optical properties		aerosol–cloud interactions		canopy nitrogen profile	atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration
			Impact of BBA	More	More	Aerosols only	No AIE	Steeper nitrogen	$+\mathrm{25}$ ppm of	$+\mathrm{50}$ ppm of
			emissions on	scattering	absorbing	affect cloud effective		profile in	atmospheric	atmospheric
			vegetation	BBA	BBA	radius (1stAIE)		the canopy	${\mathrm{CO}}_{\mathrm{2}}$	${\mathrm{CO}}_{\mathrm{2}}$
Present-day BBA emissions	over South America	Control runs, i.e.	BBAx0	BBAx0${}_{\mathrm{DIFF}\mathrm{\_}\mathrm{OP}}$	BBAx0${}_{\mathrm{ABS}\mathrm{\_}\mathrm{OP}}$	BBAx0${}_{\mathrm{1}\mathrm{stAIE}}$	BBAx0${}_{\mathrm{noAIE}}$	BBAx0${}_{\mathrm{Steep}\mathrm{\_}\mathrm{N}}$	BBAx0${}_{+\mathrm{25}\mathrm{ppm}}$	BBAx0${}_{+\mathrm{50}\mathrm{ppm}}$
		no BBA
		emission
		Half the	BBAx0.5	–	–	–	–	–	–	–
		BBA emissions
		Standard BBA	BBAx1	BBAx1${}_{\mathrm{DIFF}\mathrm{\_}\mathrm{OP}}$	BBAx1${}_{\mathrm{ABS}\mathrm{\_}\mathrm{OP}}$	BBAx1${}_{\mathrm{1}\mathrm{stAIE}}$	BBAx1${}_{\mathrm{noAIE}}$	BBAx1${}_{\mathrm{Steep}\mathrm{\_}\mathrm{N}}$	BBAx1${}_{+\mathrm{25}\mathrm{ppm}}$	BBAx1${}_{+\mathrm{50}\mathrm{ppm}}$
		emissions
		Double the	BBAx2	BBAx2${}_{\mathrm{DIFF}\mathrm{\_}\mathrm{OP}}$	BBAx2${}_{\mathrm{ABS}\mathrm{\_}\mathrm{OP}}$	BBAx2${}_{\mathrm{1}\mathrm{stAIE}}$	BBAx2${}_{\mathrm{noAIE}}$	BBAx2${}_{\mathrm{Steep}\mathrm{\_}\mathrm{N}}$	BBAx2${}_{+\mathrm{25}\mathrm{ppm}}$	BBAx2${}_{+\mathrm{50}\mathrm{ppm}}$
		BBA emissions
		Quadruple the	BBAx4	–	–	–	–	–	–	–
		BBA emissions

Aerosol optical properties

The representation of BBA in HadGEM2-ES is based on the measurements collected during the SAFARI 2000 campaign near South Africa (Abel et al., 2003; Bellouin et al., 2011). It describes the size distribution of BBA as an external mixture of two mono-modal smoke species. For the fresh smoke, a log-normal distribution with a median geometrical radius ($r$), $r=\mathrm{0.1}$ $\mathrm{\mu }\mathrm{m}$, and a geometric standard deviation ($\mathit{\sigma })$, $\mathit{\sigma }=\mathrm{1.30}$, are assumed. For aged smoke, $r=\mathrm{0.12}$ $\mathrm{\mu }\mathrm{m}$ and $\mathit{\sigma }=\mathrm{1.30}$. Fresh biomass smoke is converted to aged smoke at an exponential rate assuming an e-folding time of 6 h, which typically accounts for the ageing of the smoke plume due to condensation of chemical species (e.g. sulfate or organic compounds; Abel et al., 2003). Optical properties for the two modes are calculated a priori (i.e. offline) using Mie theory for various levels of relative humidity (RH) to account for hygroscopic growth. These optical properties – specific extinction, absorption coefficients and asymmetry parameter – are then prescribed in the HadGEM2-ES radiative transfer look-up table of optical properties.

BBA optical properties may vary significantly depending on the type of vegetation burnt, combustion regime and the meteorological conditions (Reid et al., 2005). Many observational campaigns since SAFARI 2000 have reported more absorbing BBA in other regions of the world (e.g. Johnson et al., 2008, 2016). Even at the regional scale, variation in BBA optical properties may occur. For example, aircraft observations in Brazil during SAMBBA show that flaming combustion associated with Cerrado burning in the eastern regions produces more BC and less organic aerosol, and therefore a more absorbing BBA, while smouldering forest burning in the west produces a less absorbing BBA (Johnson et al., 2016). The degree of aerosol absorption is characterised by the single scattering albedo (SSA), which is the ratio of aerosol scattering over aerosol extinction. BBAs with low SSA (e.g. $\sim \mathrm{0.80}$) absorb more solar radiation than BBAs with higher SSA (e.g. $\sim \mathrm{0.90}$). This can have implications from the vegetation perspective as a layer made of absorbing BBA would transmit less radiation to the surface than a layer made of a more scattering BBA, limiting the amount of energy available for photosynthesis. In this experiment, we investigate this aspect by varying BBA SSA by $\pm \mathrm{10}$ % by scaling the specific scattering ($K_{\mathrm{sca}}$ in m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$) and absorption ($K_{\mathrm{abs}}$ in m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$) coefficients ($K_{\mathrm{sca}}$ in m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$) directly in the look-up tables, ensuring that specific extinction remains constant. The asymmetry parameter is assumed to be unaffected. Dry BBA optical properties at 550 nm for the aged smoke are reported in Table 2.

Dry (relative humidity is 0 %) optical properties at 550 nm for the aged smoke biomass burning aerosols.

	$K_{\mathrm{ext}}$ (m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$)	$K_{\mathrm{abs}}$ (m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$)	$K_{\mathrm{sca}}$ (m${}^{\mathrm{2}}$ kg${}^{-\mathrm{1}}$)	SSA
Scattering BBA	$\mathrm{5.073}\times {\mathrm{10}}^{\mathrm{3}}$	$\mathrm{9.191}\times {\mathrm{10}}^{\mathrm{2}}$	$\mathrm{4.154}\times {\mathrm{10}}^{\mathrm{3}}$	0.99
Standard BBA	$\mathrm{5.073}\times {\mathrm{10}}^{\mathrm{3}}$	$\mathrm{4.575}\times {\mathrm{10}}^{\mathrm{2}}$	$\mathrm{4.615}\times {\mathrm{10}}^{\mathrm{3}}$	0.91
Absorbing BBA	$\mathrm{5.073}\times {\mathrm{10}}^{\mathrm{3}}$	$\mathrm{5.074}\times {\mathrm{10}}^{-\mathrm{1}}$	$\mathrm{5.072}\times {\mathrm{10}}^{\mathrm{3}}$	0.82

For this sensitivity experiment, the BBAx0, BBAx1 and BBAx2 simulations are rerun twice, once assuming a more absorbing BBA and once assuming a more scattering BBA (simulations labelled BBAx0${}_{\mathrm{DIFF}\mathrm{\_}\mathrm{OP}}$, BBAx1${}_{\mathrm{DIFF}\mathrm{\_}\mathrm{OP}}$ and BBAx2${}_{\mathrm{DIFF}\mathrm{\_}\mathrm{OP}}$ for the diffuse case and BBAx0${}_{\mathrm{ABS}\mathrm{\_}\mathrm{OP}}$, BBAx1${}_{\mathrm{ABS}\mathrm{\_}\mathrm{OP}}$ and BBAx2${}_{\mathrm{ABS}\mathrm{\_}\mathrm{OP}}$ for the absorbing case, respectively). Figure S6 in the Supplement shows how HadGEM2-ES simulates the ambient SSA of BBA (Fig. S6a) and of all aerosols (Fig. S6b) after modifying the BBA optical properties. Figure S6c shows that the amount of direct PAR is unaffected as expected because of the constraint imposed on $K_{\mathrm{ext}}$. In the higher SSA case (i.e. more diffusing BBA), the amount of diffuse PAR reaching the surface is increased, resulting in a higher amount of total PAR which contrasts with the lower SSA case.

Aerosol–cloud interactions

Clouds critically affect the amount of radiation reaching the surface (e.g. Roderick et al., 2001; Cohan et al., 2002; Pedruzo-Bagazgoitia et al., 2017). Aerosols have the potential to alter cloud properties (i.e. how they interact with radiation; Haywood and Boucher, 2000) and hence alter surface radiation. This experiment aims to address whether aerosols can affect vegetation productivity indirectly by interacting with clouds. Although aerosol–cloud interactions remain very challenging to represent in ESMs (Ghan et al., 2016; Malavelle et al., 2017), we will investigate whether the representation of these processes in the ESM used here can have a detectable impact over the region considered in this study. The BBAx0, BBAx1 and BBAx2 simulations are done twice. In the first set of simulations (labelled BBAx0${}_{\mathrm{1}\mathrm{stAIE}}$, BBAx1${}_{\mathrm{1}\mathrm{stAIE}}$ and BBAx2${}_{\mathrm{1}\mathrm{stAIE}}$), aerosols impact on precipitation efficiency is switched off (i.e. no second aerosol indirect effect, 2ndAIE, through alteration of liquid water path via auto-conversion) but can still modify cloud albedo by altering the cloud droplet effective radius (i.e. the first aerosol indirect effect, 1stAIE). In the second set of simulations (labelled BBAx0${}_{\mathrm{noAIE}}$, BBAx1${}_{\mathrm{noAIE}}$ and BBAx2${}_{\mathrm{noAIE}}$), all aerosol indirect effects are switched off. As turning off AIE reverts back CDNC to prescribed values, the BBA effect on vegetation will be calculated as a difference between simulations with the same indirect effect configuration (e.g. BBAx1${}_{\mathrm{1}\mathrm{stAIE}}$–BBAx0${}_{\mathrm{1}\mathrm{stAIE}}$).

Canopy nitrogen profile

Photosynthesis not only requires light, ${\mathrm{CO}}_{\mathrm{2}}$ and water but also nutrients that are essential in the chemistry cycles of photosynthesis. Nitrogen can be considered the most critical of those nutrients and could act as a bottleneck for plant photosynthesis (e.g. Bonan, et al., 2011; Ciais et al., 2014; Fernández-Martínez et al., 2014; Wieder et al., 2015; Houlton et al., 2015; Zaehle et al., 2015). Optimisation arguments suggest that, in order to maximise the rates of carboxylation and the rate of transport of photosynthetic products, nitrogen resources should be allocated at the top of the canopy (i.e. a steep decrease in the nitrogen profile) where light absorption is maximum (Alton and North, 2007). However, observations support a more even allocation of the nitrogen resources (i.e. a shallow decrease in the nitrogen profile throughout the canopy; Mercado et al., 2009b; Lloyd et al., 2010; Dewar et al., 2012).

Nitrogen limitation and the nitrogen cycle are not yet represented explicitly in HadGEM2-ES but will be in future versions of this Earth system model (i.e. UKESM1). Presently, nitrogen allocation at the leaf level (N${}_{\mathrm{Leaf}}$) within the canopy is represented via an exponential profile in the land surface code of HadGEM2-ES, that is $${\mathrm{N}}_{\mathrm{Leaf}}\left(L\right)={\mathrm{N}}_{\mathrm{L}\mathrm{0}}{e}^{-K_{\mathrm{N}}\mathrm{L}},$$ where $L$ is the leaf-level leaf area index, N${}_{\mathrm{L}\mathrm{0}}$ is the nitrogen concentration at canopy top (in kgN kgC${}^{-\mathrm{1}}$) and $K_{\mathrm{N}}$ is a dimensionless constant representing the steepness of the nitrogen profile. A shallow nitrogen profile ($K_{\mathrm{N}}=\mathrm{0.128}$) is the JULES default (Mercado et al., 2007) and is assumed in HadGEM2-ES for the main experiment. For this sensitivity experiment, we investigate the consequence of assuming a steeper nitrogen profile ($K_{\mathrm{N}}=\mathrm{0.5}$). Under these conditions, one might expect lesser light use efficiency under diffuse light conditions as shaded leaves become nitrogen limited (Hikosaka, 2014). We rerun the BBAx0, BBAx1 and BBAx2 simulations using the steeper nitrogen profile (labelled BBAx0${}_{\mathrm{STEEP}\mathrm{\_}\mathrm{N}}$, BBAx1${}_{\mathrm{STEEP}\mathrm{\_}\mathrm{N}}$ and BBAx2${}_{\mathrm{STEEP}\mathrm{\_}\mathrm{N}}$, respectively).

To derive a new parameter value of $K_{\mathrm{N}}$ which still provides consistent global NPP fluxes, we repeated the offline analysis described in Sect. 2.1. We used JULES to perform 1-D simulations of a tropical site with varying combinations of the $K_{\mathrm{N}}$ and N${}_{\mathrm{L}\mathrm{0}}$ parameters to derive biochemical fluxes (Fig. S4b–c). The parameter combinations were chosen such that the mean canopy carboxylation rate ($V_{\mathrm{cmax},\mathrm{25},\mathrm{C}}$) is conservative and remained at the same level as in the main experiment (i.e. about 27 $\mathrm{\mu }\mathrm{mol}$ ${\mathrm{CO}}_{\mathrm{2}}$ m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for broadleaf trees). With nitrogen allocation being represented by an exponential decay, the mean canopy $V_{\mathrm{cmax},\mathrm{25},\mathrm{C}}$ can be calculated as follows: $${\displaystyle }{V}_{\mathrm{cmax},\mathrm{25},\mathrm{C}}=n_e{N}_{\mathrm{L}\mathrm{0}}{\displaystyle \frac{\left(\mathrm{1}-e^{-K_{\mathrm{N}}\mathrm{LAI}}\right)}{K_{\mathrm{N}}\mathrm{LAI}}},$$ where LAI is the leaf area index at canopy level, $n_e$ is a constant that has values of 0.0008 and 0.0004 mol ${\mathrm{CO}}_{\mathrm{2}}$ m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ kgC (kgN)${}^{-\mathrm{1}}$ for C${}_{\mathrm{3}}$ and C${}_{\mathrm{4}}$ plants, respectively (Mercado et al., 2007).

Atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration

It is hypothesised that in a richer ${\mathrm{CO}}_{\mathrm{2}}$ world, rates of photosynthesis would increase and in addition plants could afford a reduced stomatal opening to fix the same amount of ${\mathrm{CO}}_{\mathrm{2}}$, resulting in a higher water use efficiency which should further enhance plant productivity – the so-called ${\mathrm{CO}}_{\mathrm{2}}$ fertilisation effect (e.g. Keenan et al., 2013). As stated earlier, fires do not only release aerosol particles but also ${\mathrm{CO}}_{\mathrm{2}}$, amongst other gases, which locally increases background ${\mathrm{CO}}_{\mathrm{2}}$ levels (e.g. Wittenberg et al., 1998). Additionally, it is expected that the rise in atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ will continue given current projections of anthropogenic emissions (O'Neill et al., 2016). The details of the ${\mathrm{CO}}_{\mathrm{2}}$ fertilisation effect are complex because environmental changes occur simultaneously (e.g. van der Sleen et al., 2015; Zhu et al., 2016). It would be far beyond the scope of this study to fully characterise the ${\mathrm{CO}}_{\mathrm{2}}$ fertilisation effect strength in HadGEM2-ES, but it is certainly of interest to evaluate if the effect of aerosols on vegetation through alteration of the surface PAR differs when the atmospheric background ${\mathrm{CO}}_{\mathrm{2}}$ is varied. For this experiment, the BBAx0, BBAx1 and BBAx2 simulations are done twice: once with the level of background ${\mathrm{CO}}_{\mathrm{2}}$ increased by $+\mathrm{25}$ ppm globally and once with an increase of $+\mathrm{50}$ ppm globally. Increments of $+\mathrm{25}$ and $+\mathrm{50}$ ppm should be representative of the ${\mathrm{CO}}_{\mathrm{2}}$ level expected in 12.5 and 25 years, respectively, if one assumes a 2 ppm yr${}^{-\mathrm{1}}$ increase (as supported by the mean rate of ${\mathrm{CO}}_{\mathrm{2}}$ increase measured at Mauna Loa for the period 2000–2010).

A framework to analyse the changes in fast carbon fluxes

As stated previously, aerosols can affect photosynthetic rates through different pathways (e.g. Bonan, 2008 and Fig. S7). Firstly, by altering the amount of light (the reduction in total PAR) and light quality (the change in diffuse fraction of PAR). Secondly, aerosols interact with radiation and clouds impacting the climate directly and indirectly which affects the radiative balance therefore the energy budget, forcing the coupled system to adjust to the aerosol perturbations. These adjustments (the climate feedback) can feedback into the calculations of the rate of vegetation biochemical processes – e.g. by altering the surface temperature. A simple theoretical framework can be used to discriminate a fast carbon flux, e.g. NPP, as a function of the diffuse fraction, $f_{\mathrm{d}}$, the total PAR, TotPAR, and the climate feedback, clim, such as NPP($f_{\mathrm{d}}$, TotPAR, clim). Neglecting the interdependency between the three terms enables the following decomposition: $${\displaystyle }\mathit{\delta }\mathrm{NPP}\cong {\displaystyle \frac{\partial \mathrm{NPP}}{\partial f_{\mathrm{d}}}}\mathit{\delta }{f}_{\mathrm{d}}+{\displaystyle \frac{\partial \mathrm{NPP}}{\partial \mathrm{TotPAR}}}\mathit{\delta }\mathrm{TotPAR}+{\displaystyle \frac{\partial \mathrm{NPP}}{\partial \mathrm{Clim}}}\mathit{\delta }\mathrm{Clim}.$$ To evaluate how these three terms contribute individually to the total change in NPP (the net impact), we have developed three new model diagnostics in HadGEM2-ES. For each model time step, we diagnose four surface fluxes of PAR which are the total and direct PAR, considering or excluding the aerosol radiative effects. This is achieved by calling the radiative transfer routines twice (i.e. a double call) within the same model time step, i.e. first call with the aerosol radiative effects considered and second call assuming “clean-sky” conditions where the radiative effects of aerosols are not considered (Ghan, 2013). Note that the effect of clouds on the radiative fluxes are always considered during the two calls. The next model iteration (i.e. the prognostic call) always includes the aerosol radiative effects in order to account for their impact on the atmospheric state. That means that the calculation of vegetation processes which occurs after the radiative transfer will always “see” the climate that has been modified by the aerosols. After the radiative transfer calculations, the four fluxes of PAR that were calculated are passed to the physiology routines of JULES to calculate plant productivity. Prior to calculating the biochemical fluxes, we define two values of $f_{\mathrm{d}}$ and TotPAR using the four PAR fluxes previously introduced; one that considers the effect of aerosols ($f_{\mathrm{d}}$.aer and TotPAR.aer) and one that considers clean-sky conditions ($f_{\mathrm{d}}$.clean and TotPAR.clean).

Model quantities calculated during the triple call of the physiology routines (see text).

			Aerosol effect on model variables during the triple call:
			with (.aer) and without (.clean) aerosol effect
			$f_{\mathrm{d}}$	TotPAR	clim	Biochemical flux diagnostic (e.g. NPP)	Comments
Call order of the	physiology routines	No. 1	$f_{\mathrm{d}}$.clean	TotPAR.clean	clim.aer	NPP${}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{clean},f{}_{\mathrm{d}}.\mathrm{clean}}$	NPP of vegetation
							only experiencing the
							change in climate
		No. 2	$f_{\mathrm{d}}$.clean	TotPAR.aer	clim.aer	NPP${}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{aer},f{}_{\mathrm{d}}.\mathrm{clean}}$	No. 2 minus no. 1
							is the impact of change in
							total amount of PAR
		No. 3	$f_{\mathrm{d}}$.aer	TotPAR.aer	clim.aer	NPP${}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{clean},f{}_{\mathrm{d}}.\mathrm{aer}}$	No. 3 minus no. 2
							is the impact of change in
							diffuse fraction of PAR

The physiology routines are then called three times (i.e. a triple call, see Table 3) within the same model time step. On the first call, both the reduction in total PAR and the change in diffuse fraction are ignored (i.e. the vegetation only sees the climate feedback). The biochemical fluxes calculated during this first call are saved in a specific model diagnostic (NPP${}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{clean},f_{\mathrm{d}}.\mathrm{clean}}^{\mathrm{BBAxx}}$). On the second call, the reduction in total PAR due to aerosols is then considered, but the change in diffuse fraction of PAR is not accounted for, and a new set of biochemical fluxes are saved in a specific model diagnostic (NPP${}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{aer},f_{\mathrm{d}}.\mathrm{clean}}^{\mathrm{BBAxx}}$). For the last prognostic call, both aerosol effects on reduction in total PAR and the change in diffuse fraction are taken into account in the calculation of the biochemical fluxes and saved in a specific model diagnostic (NPP${}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{aer},f_{\mathrm{d}}.\mathrm{aer}}^{\mathrm{BBAxx}}$).

With these new diagnostics available, we are able to isolate the impacts of change in diffuse fraction, reduction in total PAR and climate feedback by comparing model simulations which include or exclude the BBA emissions. For instance, the effect of BBA in the BBAx1 simulation (i.e. the standard emissions scenario) can be expressed as follows: $$\begin{array}{cc}& {\displaystyle }\mathrm{\Delta }{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{net}\mathrm{impact}}^{\mathrm{BBAx}\mathrm{1}}& {\displaystyle }={\stackrel{\mathrm{‾}}{\mathrm{NPP}}}^{\mathrm{BBAx}\mathrm{1}}-{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}^{\mathrm{BBAx}\mathrm{0}}{\displaystyle }& {\displaystyle }\cong \mathrm{\Delta }{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{f_{\mathrm{d}}}^{\mathrm{BBAx}\mathrm{1}}+\mathrm{\Delta }{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{TotPAR}}^{\mathrm{BBAx}\mathrm{1}}+\mathrm{\Delta }{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{clim}}^{\mathrm{BBAx}\mathrm{1}}\end{array}$$ with $$\begin{array}{cc}{\displaystyle }\mathrm{\Delta }{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{f_{\mathrm{d}}}^{\mathrm{BBAx}\mathrm{1}}=& {\displaystyle }\left({\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{aer},f_{\mathrm{d}}.\mathrm{aer}}^{\mathrm{BBAx}\mathrm{1}}\right.\\ {\displaystyle }& {\displaystyle }\left.-{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{aer},f_{\mathrm{d}}.\mathrm{clean}}^{\mathrm{BBAx}\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }-\left({\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{aer},f_{\mathrm{d}}.\mathrm{aer}}^{\mathrm{BBAx}\mathrm{0}}\right.\\ & {\displaystyle }& {\displaystyle }-\left.{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{aer},f_{\mathrm{d}}.\mathrm{clean}}^{\mathrm{BBAx}\mathrm{0}}\right),{\displaystyle }\mathrm{\Delta }{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{TotPAR}}^{\mathrm{BBAx}\mathrm{1}}=& {\displaystyle }\left({\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{aer},f_{\mathrm{d}}.\mathrm{clean}}^{\mathrm{BBAx}\mathrm{1}}\right.\\ {\displaystyle }& {\displaystyle }-\left.{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{clean},f_{\mathrm{d}}.\mathrm{clean}}^{\mathrm{BBAx}\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }-\left({\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{aer},f_{\mathrm{d}}.\mathrm{clean}}^{\mathrm{BBAx}\mathrm{0}}\right.\\ & {\displaystyle }& {\displaystyle }\left.-{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{clean},f_{\mathrm{d}}.\mathrm{clean}}^{\mathrm{BBAx}\mathrm{0}}\right),{\displaystyle }\mathrm{\Delta }{\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{clim}}^{\mathrm{BBAx}\mathrm{1}}=& {\displaystyle }\left({\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{Clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{clean},f_{\mathrm{d}}.\mathrm{clean}}^{\mathrm{BBAx}\mathrm{1}}\right)\\ & {\displaystyle }& {\displaystyle }-\left({\stackrel{\mathrm{‾}}{\mathrm{NPP}}}_{\mathrm{Clim}.\mathrm{aer},\mathrm{TotPAR}.\mathrm{clean},{\mathrm{f}}_{\mathrm{d}}.\mathrm{clean}}^{\mathrm{BBAx}\mathrm{0}}\right),\end{array}$$ where overbars denote quantities averaged over a time period long enough for vegetation fast responses to adjust to the aerosol effects.

Global annual estimates of gross primary productivity (GPP, a, c, e) and net primary productivity (NPP, b, d, f). Observationally based estimates from FLUXCOM MTE analysis (a), MODIS MOD17A2 (b) and HadGEM2-ES (c, d). Zonal means are shown in panels (e) and (f). The circles on the NPP maps (b, d) represent in situ estimates from the EMDI project.

[Figure omitted. See PDF]

Observations used in model evaluation

We evaluate global fields of simulated GPP and NPP using GPP fields derived by the FLUXCOM project (Tramontana et al., 2016; Jung et al., 2017a) and the global annual mean NPP retrievals based on the MODIS MOD17A2 product (Running et al., 1994) (Fig. 1a, b). The GPP from FLUXCOM is derived from a model that has been trained on observational data, so we will refer to this estimate as a “reconstructed” GPP. In addition, in situ estimates of NPP from the EMDI project (http://gaim.unh.edu/Structure/Intercomparison/EMDI/, last access: 25 January 2019) are also presented in the form of overlaid circles depicted in Fig. 1b. Note, simulated values of HadGEM2-ES GPP and NPP used in the comparison with observational data are sampled where the corresponding observationally based dataset contains non-missing data.

The simulated aerosol loading is evaluated against the record of aerosol optical thicknesses (AOTs) retrieved from the MODIS instrument measurements on board of the Terra satellite. The dataset used corresponds to the Level-3 MODIS Atmosphere Monthly Global Product collection 6.1 (at 1-degree resolution) that was derived from the MYD06_L2 products for the period extending between 2001 and 2016.

Additional evaluation of the model skill against observations is provided in the Supplement (Fig. S8). This includes comparisons of the modelled solar fluxes at the surface against the SSF1deg Terra Edition 2.8 product based on the CERES radiation data, and comparisons of the modelled surface precipitation against the GPCP version 2.3 product.

Multi-annual mean for the June–July–August season (JAS) of the aerosol optical thickness (AOT) at 550 nm (a, b) and the seasonal cycle (c, d) of the AOT calculated over the domain highlighted in red for the MODIS Terra retrieval (a, c) and the HadGEM2-ES model (b, d). The MODIS seasonal cycle (c) shows the multi-year (2001–2016) mean with the black line, and the individual years are overlaid with red dashed lines. The seasonal cycle for HadGEM2-ES (d) shows the 30-year mean for the five experiments with varying biomass burning emissions (see text, Sect. 2.2).

[Figure omitted. See PDF]

Results

Evaluation

Carbon exchange

Global annual mean GPP and NPP as simulated by HadGEM2-ES with the new representation of canopy light interception are shown in Fig. 1c, d. The global GPP modelled by HadGEM2-ES is $+\mathrm{115}$ PgC yr${}^{-\mathrm{1}}$ in the updated version of HadGEM2-ES and smaller than the estimate of $+\mathrm{129}$ PgC yr${}^{-\mathrm{1}}$ from the FLUXCOM dataset (Fig. 1a) but closer to the reference of $+\mathrm{118}$ PgC yr${}^{-\mathrm{1}}$ cited by Shao et al. (2013). The standard configuration of HadGEM2-ES that participated in CMIP5 had a global GPP of the order of $+$140 PgC yr${}^{-\mathrm{1}}$ for present-day conditions (Shao et al., 2013). The underestimation of the GPP in the updated HadGEM2-ES configuration is comparable in magnitude to the overestimation of the GPP in the HadGEM2-ES configuration. However, the ratio of NPP over GPP (not shown) in the updated version of HadGEM2 is more consistent with observationally based ratio estimates (e.g. Luyssaert et al., 2007). Despite the inherent uncertainties in the two reference estimates of the global GPP (i.e. between $+$118 and 129 TgC yr${}^{-\mathrm{1}}$), it suggests that the updated version of HadGEM2-ES is able to provide a more consistent global GPP estimate. Over the central Amazon domain, which is represented by the region encapsulated in the red box in Fig. 2a., the HadGEM2-ES average GPP in August (respectively September) is $\mathrm{2750}\pm \mathrm{250}$ gC m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ (respectively $\mathrm{2600}\pm \mathrm{200}$ gC m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ for September) compared to $\mathrm{2250}\pm \mathrm{125}$ gC m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ (respectively $\mathrm{2500}\pm \mathrm{180}$ gC m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for September) for FLUXCOM.

The global NPP modelled by HadGEM2-ES is $+\mathrm{54}$ PgC yr${}^{-\mathrm{1}}$ (Fig. 1d) and in good agreement with the satellite-based estimate of $+\mathrm{50}$ PgC yr${}^{-\mathrm{1}}$ (Fig. 1b) and the “best guess” value of $+\mathrm{56}$ PgC yr${}^{-\mathrm{1}}$ reported by Shao et al. (2013). The updated configuration of HadGEM2-ES performs well in mid and high latitudes, particularly against EMDI data (Fig. 1d), but biases still remain in the tropics (Fig. 1f), particularly over South America in areas dominated by C3 grass (Fig. S5). Despite obvious overestimation by HadGEM2-ES of the NPP on annual mean over South America when compared to MODIS MOD17A2 (Fig. 1b, d), the fluxes are well captured during the peak of the fire season over the central Amazon. The average GPP from HadGEM2-ES in August (respectively September) is $\mathrm{1080}\pm \mathrm{140}$ gC m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ (respectively $\mathrm{975}\pm \mathrm{100}$ gC m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ for September) compared to $\mathrm{990}\pm \mathrm{550}$ gC m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ (respectively $\mathrm{1025}\pm \mathrm{590}$ gC m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for September) for MODIS MOD17A2.

Biomass burning aerosols

Biomass burning is highly variable from year to year. This can be readily observed by monitoring the AOT, a proxy for the amount of aerosol particles present in the atmosphere. Figure 2a shows the average AOT retrieved at 550 nm for the months July–August–September (JAS) between 2001 and 2016 by the MODIS instrument on board of the Terra satellite. Although most of man-made fires occur in the so-called arc of deforestation on the edge of the rainforest (Cochrane, 2003), the hot spot of high AOT (> 0.6) is actually observed over the Rondônia state (Brazil) near the Bolivian border. This hotspot can be explained by (i) the action of the large-scale atmospheric circulation that recirculates aerosols over South America and (ii) the contribution of natural fires that occur concomitantly with fires of anthropogenic origin. Figure 2c provides more detail on the AOT variability by showing the seasonal cycle calculated over the central Amazon (i.e. the region encapsulated in the red box shown in Fig. 2a using the multi-year data record from MODIS). Despite year-to-year variability, AOT is found to peak in September over this region that is, at the expected peak of the fire season, supporting that BBAs are the dominant component of the total aerosol loading during that period.

The AOT modelled by HadGEM2-ES in the simulation that assumes standard BBA emission (i.e. the BBAx1 simulation) is in overall good agreement with the MODIS observations for the JAS period (Fig. 2a, b; Johnson et al., 2016). However, the AOT at the peak of the fire season (i.e. in September) is underestimated (Fig. 2d). In contrast, the modelled AOT for September in the BBAx2 simulation is in better agreement with the satellite retrievals. We will therefore consider in the remainder of this paper that the combination of BBAx1 and BBAx2 scenarios are representative of present-day levels of BBA and will use them to discuss the effects of BBA on the rainforest productivity. There is huge variation in the inter-annual variation in the magnitude of the AOT (Fig. 2c), which justifies the upper bound for our simulation scenarios; the simulations BBAx0.5 and BBAx4 will be considered representative of emissions for years with low and high fire activity, respectively (Fig. 2c). These simulations will provide a lower and upper estimate, respectively, of the BBA impact on vegetation.

Modelled seasonal cycle from HadGEM2-ES for the total PAR (a, b), the diffuse PAR (c, d) and fraction of radiation that is diffuse (e, f) for the five BBA emission experiments. Absolute values (a, c, e) and relative anomalies (b, d, f) $w/r$ to experiment BBAx0 (i.e. no biomass burning aerosols) are shown. Transparent coloured areas in panels (a, c, e) correspond to $\pm \mathrm{1}$ standard deviation. Dashed lines are the multi-year annual means.

[Figure omitted. See PDF]

Surface radiation

Figure 3 illustrates the impact of BBA on the radiative fluxes in the HadGEM2-ES simulations. The seasonal cycle of the total PAR (TotPAR) shows a strong decrease during the whole dry season with the strongest reduction occurring in August and September. The reduction in TotPAR is in the range of $-\mathrm{18.0}$ to $-\mathrm{7.5}$ W m${}^{-\mathrm{2}}$ (i.e. $-\mathrm{14.0}$ to $-\mathrm{5.5}$ %) in the BBAx1 and BBAx2 experiments, respectively (Fig. 3a, b). For the most extreme emission scenarios (BBAx4), the reduction in TotPAR is as high as $-\mathrm{30}$ W m${}^{-\mathrm{2}}$, or $-\mathrm{25}$ %, in August. Conversely, the diffuse component of PAR (DiffPAR) increases with aerosols as expected from the theory of light scattering (Fig. 3c, d). The diffuse PAR reaching the top of the canopy is increased by approximately $+\mathrm{6.0}$ to $+\mathrm{12.0}$ W m${}^{-\mathrm{2}}$ (i.e. approximately $+\mathrm{14.0}$ to $+\mathrm{31.0}$ %) during August and September in the BBAx1 and BBAx2 simulations (Fig. 3c, d). Overall this leads to an increase in the diffuse fraction of PAR (i.e. $f_{\mathrm{d}}$) of $+\mathrm{20.0}$ to $+\mathrm{55.0}$ % (Fig. 3e, f).

Showing the total PAR (TOTPAR, a), the diffuse PAR (DIFFPAR, b) and the fraction of PAR that is diffused (DIFF_FRAC, c) reaching the surface versus the total aerosol optical thickness (AOT) at 550 nm and the net primary productivity (NPP, d) against the fraction of PAR that is diffused. Circles represent the binned data from the HadGEM2-ES simulations, while plain lines are the corresponding second-order polynomial fits. Prior to binning, data were first collected at all grid cells in the Amazon region (i.e. the red box region in Fig. 2) for all five BBA emission experiments. We then aggregate all grid cells into 30 AOT bins ranging from 0 to 3 at an interval of 0.1. In each bin, we calculate average AOT and corresponding TOTPAR, DIFFPAR and DIFF_FRAC (Fig. 4a, b and c, respectively; we calculate average DIFF_FRAC and corresponding NPP in Fig. 4d).

[Figure omitted. See PDF]

An alternative representation of the impact of BBA on the radiative fluxes is depicted in Fig. 4 for August and September. Here, the composite plot is constructed using the four simulations that include BBA emissions to calculate the TotPAR (Fig. 4a), DiffPAR (Fig. 4b) and $f_{\mathrm{d}}$ (Fig. 4c) at the surface as a function of the total AOT (i.e. BBA $+$ background aerosols). The composite was constructed by first averaging each simulation over time to create climatologies for the specific months, and then all pixels contained in the domain of analysis were sampled to construct the scatterplots of the desired quantities. It is important to note that radiative quantities were sampled for the full sky grid box and that no conditional sampling was applied a priori; therefore cloud effects are implicitly accounted for in these statistics. Subsequently, further averaging of the data into 30 bins of AOT (respectively $f_{\mathrm{d}}$ for Fig. 4d) was applied to smooth the signal. Figure 4a shows the expected monotonic decrease in TotPAR with AOT. Concomitantly, the DiffPAR (Fig. 4b) increases with AOT up to values of around 1.75 and decreases for higher AOTs. This illustrates that increasing AOT could only increase the amount of diffuse light reaching the surface up to a point; above this point, the effect of the attenuation of TotPAR dominates. This AOT threshold around 1.75 maximises the amount of diffuse radiation reaching the canopy top. However, as it will be detailed in following sections, this threshold does not correspond to the maximum effect of aerosols on vegetation productivity.

The net impact of BBA on forest productivity

Figure 4d represents NPP as a function of $f_{\mathrm{d}}$ for the months of August and September in the same way as the surface radiative fluxes against AOT are depicted (Fig. 4a–c). This shows that NPP is likely to reach an optimum when $f_{\mathrm{d}}$ approximately equals 52 %–56 %. The existence of an optimum $f_{\mathrm{d}}$ that would maximise carbon sequestration is consistent with findings reported in past modelling studies (e.g. Knohl and Baldocchi, 2008; Mercado et al., 2009a; Pedruzo-Bagazgoitia et al., 2017; Yue and Unger, 2017). Such an optimum, however, depends strongly on factors such as the vegetation canopy architecture environmental conditions, solar zenith angle or the optical properties of the scattering particles. The fact that an optimum diffuse fraction emerges is consistent with our understanding of the DFE mechanism. When $f_{\mathrm{d}}$ is lower than the optimum, an increase in the amount of diffuse radiation increases carbon assimilation because a larger area of shaded leaves become photosynthetically active. For $f_{\mathrm{d}}$ beyond the optimum, the effect of the attenuation of TotPAR dominates, and sunlit leaves are no longer light saturated, resulting in an overall decrease in biochemical fluxes at the canopy level with further increase in $f_{\mathrm{d}}$.

Figure 4c could be used to infer an AOT for which $f_{\mathrm{d}}$ is getting close to the optimum value of 0.55 (Fig. 4d). This would approximately occur at an AOT of $\sim \mathrm{0.9}$–1 (Fig. 4c). However, we do not observe that the highest NPP enhancement occurs around these values of AOT in our simulations (see Sect. 3.3). This can be understood as a consequence of equifinality because both the effects of clouds and the effects of aerosols on radiation occur concomitantly. There are then many possible combinations of cloud and aerosol scenarios that could create optimum conditions maximising the DFE. It would be possible to disentangle the effect of BBA from the effect of clouds on carbon sequestration by either screening out cloudy scenes or diagnosing the biochemical fluxes in the clear-sky portion of the model grid boxes, providing a mean to quantify the maximum potential impact of BBA on carbon sequestration. This approach was used by Moreira et al. (2017) to conclude that BBA could increase the GPP of the Amazon forest by up to 27 %. While this study is insightful, our aims here are different as we seek to understand the impact of BBA while considering the system-wide behaviour that is including the effects of both aerosols and clouds. This alternative approach was used by Yue and Unger (2017) to analyse aerosol impacts on vegetation over China and show that clouds are a dominant feature, controlling the diffuse fraction of radiation which modulates the diffuse fertilisation effect from aerosols (Yue and Unger, 2017, Fig. 5 therein). In Sect. 3.3, we will show that similar conclusions could be drawn over South America.

NPP anomalies (relative to the experiment BBAx0) for the 30-year mean for the four varying BBA emissions (see text, Sect. 2.2) during the August (a, c, e, g) and September (b, d, f, h) months. Mean fluxes (labelled AVG) and accumulation (labelled TOT) are calculated over the domain delimited by the pink borders. Hatched areas represent the regions where changes are significant at the 95 % confidence level. Green contours show the 550 nm AOT anomalies.

[Figure omitted. See PDF]

Despite cloudiness affecting how much aerosols can interact with radiation, we notice that NPP is enhanced in the central part of the Amazon when BBA emissions are increased (Fig. 5). The most statistically significant enhancement of the NPP, which is depicted by the stippling in Fig. 5, occurs during August, in phase with the period when the radiative impacts of BBA are the most pronounced in the model simulations (Fig. 3, Sect. 3.1.3). Although the simulated AOTs are of similar magnitude during September, NPP enhancement is not as robust as in August (i.e. there is a less statistically significant signal in the NPP anomalies). This can be partially explained by the fact that plant productivity simulated by HadGEM2-ES reaches a minimum in September (Fig. S8a, b). As a result, the vegetation is less active in September and the potential impact of the BBA perturbation is reduced.

Mean seasonal cycle of NPP (a), relative changes (b) and absolute changes (c) for the five BBA emission scenarios (see text, Sect. 2.2) averaged over the Amazon Basin. Differences are calculated with regard to experiment BBAx0. The short-dash curves in panel (c) correspond to the accumulated anomalies (right $y$ axis).

[Figure omitted. See PDF]

Overall, based on the BBAx1 and BBAx2 simulations, we estimate that BBAs increase NPP by about $+\mathrm{80}$ to $+\mathrm{105}$ TgC yr${}^{-\mathrm{1}}$, or 1.9 % to 2.7 %, (Fig. 6b, c) over the domain of analysis. This estimate of the enhancement in carbon uptake is remarkably similar to the estimate provided by Rap et al. (2015), who found that Amazonian fires increase NPP by 1.4 %–2.8 % corresponding to an increase of $+$78 to $+$156 TgC yr${}^{-\mathrm{1}}$. This is encouraging as the authors used the stand-alone version of JULES (i.e. the land surface component in the HadGEM family of models). However, as it will be discussed in Sect. 3.3 and Sect. 4.2, we attribute the enhancement in carbon sequestration to different mechanisms. The Rap et al. (2015) study used a combination of offline models which do not account for climatic adjustment to the aerosol radiative perturbation. This supports the fact that the increase in modelled NPP results from DFE in their simulations. Conversely, we will show (Sect. 3.3) that DFE is negligible over the region considered in our model simulations, but the overall aerosol impacts on vegetation remains significant thanks to the contribution of climate feedbacks that are experienced by the vegetation.

Similar to Fig. 6c, showing the variation in NPP due solely to (a) change in diffuse fraction, (b) reduction in total PAR and (c) the climate feedback.

[Figure omitted. See PDF]

Disentangling the impact of radiation changes from those of climate adjustments

We have quantified the net impact of BBA on NPP in the previous section. Following the framework described in Sect. 2.4, we now separately address the individual contribution from the change in diffuse fraction, $f_{\mathrm{d}}$, the reduction in total PAR, TotPAR, and the climate feedbacks to the BBA net impact on vegetation productivity. Figure 7 shows the seasonal cycle of NPP anomalies averaged over the domain of analysis (left axis) and the corresponding accumulated anomalies (right axis) for the four simulations with varying BBA emissions. The increase in NPP due to the change in diffuse fraction is unambiguous (Fig. 7a), corresponding to an enhancement in plant net carbon uptake of $+\mathrm{65}$ to $+\mathrm{110}$ TgC yr${}^{-\mathrm{1}}$ in the BBAx1 and BBAx2 simulations, respectively. As expected, the reduction in total PAR has the opposite effect and systematically decreases NPP (Fig. 7b) with increasing negative NPP anomalies. This corresponds to a reduction in plant net carbon uptake of $-\mathrm{52}$ to $-\mathrm{105}$ TgC yr${}^{-\mathrm{1}}$ in the BBAx1 and BBAx2 simulations, respectively. The combination of the change in diffuse fraction and the reduction in total PAR effects represents the DFE. We estimate that the DFE from BBA increases the vegetation NPP by $+\mathrm{13}$ and $+\mathrm{5}$ TgC yr${}^{-\mathrm{1}}$ in the BBAx1 and BBAx2 simulations, respectively.

The impact of BBA on NPP via the DFE is in stark contrast with the increase in forest productivity, which we have discussed in the previous Sect. 3.2 describing the net impact of BBA ($+$80 to $+$105 TgC yr${}^{-\mathrm{1}}$ for the BBAx1 and BBAx2 simulations respectively). This would indicate that in our simulations the net impact of BBA on forest productivity is not mostly due to the DFE. Figure 7c shows that the climate feedback term is actually the dominant contribution and systematically increases NPP, with an enhancement of $+\mathrm{67}$ to $+\mathrm{100}$ TgC yr${}^{-\mathrm{1}}$ in the BBAx1 and BBAx2 simulations, respectively.

It is worth mentioning that the maximum impact of the change in diffuse fraction occurs during August in the BBAx4 simulation, which increases the NPP by $+\mathrm{41}$ TgC m${}^{-\mathrm{1}}$. The corresponding impact of the reduction in total PAR decreases NPP by $-\mathrm{66}$ TgC m${}^{-\mathrm{1}}$. This illustrates that for a year with intense burning, the system actually seems to shift past the point where the balance between the total and the diffuse PAR does not increase the efficiency of photosynthesis anymore (i.e. BBA DFE leads to reduction of $-\mathrm{42}$ TgC yr${}^{-\mathrm{1}}$ on an annual basis for the BBAx4 scenario). Interestingly, in this simulation, despite the negative impact on NPP from DFE, we note that the impact of climate feedback is much larger ($+\mathrm{194}$ TgC yr${}^{-\mathrm{1}}$), resulting in the net impact of BBA on the vegetation to be overall positive ($+\mathrm{151}$ TgC yr${}^{-\mathrm{1}}$).

Shown on panel (a) the relative changes in NPP ($\mathrm{\Delta }$NPP${}_{\mathrm{Net}}$, in grey), the relative changes in NPP due to the change in diffuse fraction ($\mathrm{\Delta }$NPP${}_{\mathrm{Frac}\mathrm{Diff}}$, in blue), reduction in total PAR ($\mathrm{\Delta }$NPP${}_{\mathrm{TOTPAR}}$, in red), the sum of change in diffuse fraction and reduction in total PAR ($\mathrm{\Delta }$NPP${}_{f{}_{\mathrm{d}}+\mathrm{TOTPAR}}$, in green, i.e. the DFE), and the climate feedback ($\mathrm{\Delta }$NPP${}_{\mathrm{Adjust}}$, in yellow) against the anomalies in the AOT at 550 nm for the month of August. Shown on panel (b) the relative changes in NPP due to change in diffuse fraction and reduction in total PAR ($\mathrm{\Delta }$NPP${}_{f{}_{\mathrm{d}}+\mathrm{TOTPAR}=\mathrm{DFE}}$) – i.e. the changes in NPP only due to change in surface radiation, the DFE, for August (green) and September (blue) as a function of the total AOT at 550 nm. Note this is shown against the total AOT. The dashed lines on panel (b) highlight the AOT thresholds where the DFE switches from a positive to a negative impact.

[Figure omitted. See PDF]

To compare the relative contribution of the DFE (i.e. change in diffuse fraction plus reduction in total PAR) and the climate feedbacks on vegetation NPP as the atmospheric aerosol content ramps up, Fig. 8a depicts the relative change in NPP (%) as a function of AOT for the month of August. This NPP change is further decomposed into individual contributions from the change in diffuse fraction (blue solid line), the reduction in total PAR (red solid line), the DFE (green solid line), the climate feedback (yellow solid line) and the net impact (black solid line). The resulting attribution plot shown in Fig. 8a was constructed in the same way as Fig. 4 (see Sect. 3.1), i.e. by first averaging each simulation over time, then sampling the NPP changes associated with each of the three terms in all the model grid boxes from the domain of analysis and finally aggregating the sampled quantities into 30 bins of AOT.

Overall, it is clear from Fig. 8a that BBAs enhance NPP across the entire range of AOT considered here (i.e. the net impact of BBA is strictly positive) which is consistent with the geographic distribution of anomalies displayed in Fig. 5. The impact of the change in diffuse fraction and the reduction in total PAR, respectively, consistently increases and decreases vegetation NPP, respectively, as discussed in the previous paragraph. However, the impact of DFE from the BBA (represented by the green solid line in Fig. 8a), changes its sign around AOT of $\sim \mathrm{0.75}$. At lower AOTs DFE from BBA contributes to an increase in NPP, whereas at higher AOTs it has the opposite effect. To help visualise the transition in the DFE regime, we have replotted the NPP changes due to the DFE contribution only in Fig. 8b. Here, the changes are represented for August and September and are shown against the total AOT (BBA $+$ background aerosols). It is interesting to note that the AOT threshold occurs at a smaller value in September (0.62) than in August (0.89). This suggests that the state of the climate have implication for the strength of the DFE from aerosols (e.g. via the amount of cloudiness).

As discussed in Sect. 3.1, changes in NPP due to DFE from BBA alone are calculated under all sky conditions which also account for cloud radiative effects. A plausible explanation for the observed reduction in the range of AOT creating a positive DFE would be that cloudiness increases over the analysed model domain between August and September (see Fig. S10) as the regional climate progresses towards the wet season. This is supported by the increase in $f_{\mathrm{d}}$ between August and September in the simulation that excludes BBA (i.e. black solid line in Fig. 3c). These results are consistent with those of Yue and Unger (2017) who discussed how the impact of anthropogenic aerosol DFE over China vary depending on the cloud cover which allows for smaller or larger perturbations in the radiative balance for the same atmospheric aerosol loading.

Shown in panel (a) is a box-and-whisker plot of the net primary productivity monthly means for August averaged over the central Amazon. Results are shown for the main experiment (see text, Sect. 2.2) and the four additional sensitivity experiments (see text, Sect. 2.3). Individual members of the 30-year run are represented by the green dashes. Black dots correspond to the ensemble mean. Dashed white lines are the 25th, 50th and 75th percentiles. Shown in panel (b) are the changes in NPP in each sensitivity experiments, calculated relative to their respective baseline simulation (e.g. X1${}_{+\mathrm{25}\mathrm{ppm}}$ – X0${}_{+\mathrm{25}\mathrm{ppm}}$ is the differences between the BBAx1 and BBAx0 simulations with $+\mathrm{25}$ ppm increase in ${\mathrm{CO}}_{\mathrm{2}}$ concentration).

[Figure omitted. See PDF]

Sensitivity experiments

Here, we present the results from the four additional sensitivity experiments described in Sect. 2.3. These experiments were designed to further elucidate the role BBA play in vegetation productivity while changing some of the underlying assumptions in the previous experiments which relate to (i) aerosol optical properties, (ii) aerosol–cloud interactions, (iii) the canopy nitrogen profile and (iv) atmospheric carbon dioxide concentration. Figure 9a shows a box-and-whisker plot of NPP averaged over the central Amazonia during August for all BBAx0, BBAx1 and BBAx2 simulations from the main experiment (those analysed in the Sect. 3) and from the four additional sensitivity experiments. The mean changes in NPP due to biomass burning aerosols are shown in Fig. 9b.

The results can be summarised as follows.

• Aerosol optical properties (experiments DIFF_OP and ABS_OP). The optical properties of BBA were altered in order to make the biomass burning aerosols more (DIFF_OP) or less (ABS_OP) scattering by modifying the BBA SSA (Fig. S6a, b). The mass specific extinction is invariant (see Sect. 2.3.3), which implies that for the same AOT, the direct radiation reaching the surface is also independent of the aerosol scattering and absorbing efficiency assumptions (Fig. S6c). More scattering or absorbing BBAs, respectively, increase or decrease the diffuse fraction of solar radiation reaching the surface (Fig. S6c). As a result, scattering BBAs should produce a stronger DFE and absorbing BBAs should analogously produce a weaker DFE. However, we do not observe a significant change in the modelled BBA impact on vegetation productivity for the varying BBA scattering and absorbing assumptions (Fig. 9b). In the standard simulations, the net change in NPP due to BBA is $+\mathrm{28.4}$ to 38.6 TgC month${}^{-\mathrm{1}}$ in August. For the DIFF_OP simulation (respectively ABS_OP) the net change in NPP is $+\mathrm{32.1}$ to 36.2 TgC month${}^{-\mathrm{1}}$ (respectively $+\mathrm{17.9}$ to 18.2 TgC month${}^{-\mathrm{1}}$). For September (not shown), we actually found that the ABS_OP simulation had the largest increase in NPP, which is not consistent with our assumption. In summary, the effect of BBA optical properties on NPP changes is within the noise and considered negligible. This can be explained in the light of the results discussed in Sect. 3.3, where we showed that the DFE from present-day BBA is small ($\sim +\mathrm{5}$ TgC month${}^{-\mathrm{1}}$ in August in BBAX1) for this model in this region of the world. Therefore, altering the ratio of diffuse fractions reaching the ground via the aerosol optical properties, that is modulating the magnitude of the DFE, does not have a measurable effect on vegetation productivity.

• Aerosol–cloud interactions (experiments 1stAIE and NoAIE). We have emphasised the potential role of clouds in Sect. 3.3. One could expect that increasing aerosol emissions which provide the necessary CCNs will increase cloud droplet numbers and reduce their sizes. The reduction in droplet size leads to cloud brightening (1stAIE) and possibly cloud amount (2ndAIE), which could eventually alter the surface radiation balance. We note that the impact of BBA on NPP is of similar magnitude in the main experiment and in the experiments without aerosol–cloud interactions (Fig. 9b) – i.e. neglecting ACI does not change the impact of BBA on vegetation productivity over the region considered. A possible explanation can be found in the type of the clouds that predominates in this region. We note that most of the precipitation in HadGEM2-ES stems from convective clouds. Aerosols are only coupled to the large-scale precipitation scheme in HadGEM2-ES (i.e. aerosols can only alter the properties of stratiform clouds). The absence of any impact from ACI over this region is then to be expected. Whether or not ACI can affect vegetation productivity remains a research topic for future studies, and these should focus on regions where aerosols and clouds are likely to interact as a consequence of the cloud representation in the models (e.g. Chameides et al., 1999). Alternatively, the ACI effects in the cloud representation should be revisited and improved in the models (Malavelle et al., 2017).

• Canopy nitrogen profile (experiment STEEP_N). We modified the shape of the nitrogen profile for the modelled canopy to represent a steeper decrease in leaf nitrogen content (Sect. 2.3.4). The available nitrogen to leaves decreases from the canopy top downwards. This change in leaf nitrogen allocation means that sunlit leaves have access to more resources, whereas shaded leaves tend to be more nitrogen limited (Hikosaka et al., 2014). Despite this modification in nitrogen availability, we do not observe a significant change in the modelled BBA impact on vegetation productivity. The reasons for this absence of sensitivity to nitrogen availability are similar as in the experiments testing the role of aerosol optical properties; i.e. the DFE from BBA is already too small to have a discernible impact, and reducing the allocated nitrogen in the shaded portion of the canopy only reduces its impact more.

• Atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration (experiments 25 and +50 ppm). While increasing atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration leads to an unambiguous increase in NPP (Fig. 9a), the BBA impact is of similar magnitude as in the main experiments (Fig. 9b). It may appear that the impact of BBA is somehow reduced in the $+\mathrm{25}$ ppm case compared to the main experiment and the $+\mathrm{50}$ ppm experiment. However, the level of model internal variability in NPP is too pronounced (Fig. 9a) to draw robust conclusions on the impact of a variation in ${\mathrm{CO}}_{\mathrm{2}}$ on the BBA-induced DFE. Note that the atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration increased globally. It was also allowed to affect the radiative balance, resulting in a warming climate in these two experiments. Potentially, this could increase the model's internal variability further. If one were to repeat these experiments, only the leaf-internal ${\mathrm{CO}}_{\mathrm{2}}$ concentration should be increased to avoid additional statistical noise produced in the warming climate.

Concluding remarks

From our model experiments we concluded that the diffuse PAR fertilisation effect from biomass burning aerosols in HadGEM2-ES (Sect. 3.3) is comparatively modest, amounting to between $+\mathrm{13}$ and $+\mathrm{5}$ TgC yr${}^{-\mathrm{1}}$ based on the result from the simulations BBAx1 and BBAx2. This may seem at odds with the $+\mathrm{78}$ to $+\mathrm{156}$ TgC yr${}^{-\mathrm{1}}$ estimate (assuming standard BBA emissions and 3 times the standard BBA emissions, respectively) reported by Rap et al. (2015), who used the JULES land surface model in an offline framework specifically designed to assess the DFE of biomass burning aerosols. Some differences between the two studies that could explain the apparent differences are obvious, such as for instance the fact that we are not reporting estimates for the BBA impact over the same area (i.e. our domain is smaller) or that we did not use the same aerosol properties or emission inventories. We recalculated the impact from biomass burning aerosols in our simulations over a larger domain that approximately matches the area considered by Rap et al. (2015). In this situation, we found that the net increase in NPP is about $+\mathrm{145}$ to $+\mathrm{148}$ TgC yr${}^{-\mathrm{1}}$ for the BBAx1 and BBAx2, respectively, of which only $+\mathrm{15}$ to $+\mathrm{5}$ TgC yr${}^{-\mathrm{1}}$ are attributable to the DFE. This confirms that the magnitude of the DFE from BBA effect is small, increasing plant productivity in our simulations over the Amazon forest.

Biases in the cloud amount, which are inherent of coarse model parameterisations, may affect the surface radiation and impact the magnitude of the DFE from biomass burning aerosols (and indeed all aerosols). Those uncertainties can be partially contained using an offline framework where the state of the model can be forced closer to the distribution of input observations. However, in this approach, internal consistency is lost by not allowing variability within non-linear relationships (e.g. how cloudiness is changed due to aerosol–radiation interactions, how plant dark respiration is changed due to the surface cooling). This then poses a problem and a risk of overestimating the response of a component (e.g. vegetation productivity) to a perturbation such as those introduced by aerosols. By including more complexity in a coupled framework as in the present study, we believe that our estimate of the DFE is more consistent, albeit low due to possible uncertainties or biases, and we would argue that earlier estimates of the DFE from BBA in this region (Rap et al., 2015) are probably on the high end. Nonetheless, despite showing that the DFE from BBA is not an efficient mechanism in our simulations over this region, we have demonstrated a pathway where BBA can significantly influence vegetation productivity. We assessed this pathway by calculating the term representing the biomass burning aerosol climate impact on vegetation, which represents the rapid adjustments of land surface climate to aerosol radiation perturbation. We estimated this term to be about $+\mathrm{67}$ to $+\mathrm{100}$ TgC yr${}^{-\mathrm{1}}$ over the domain analysed in this study in the BBAx1 and BBAx2 simulations, respectively. This is a novel contribution which could not be accounted for in an offline modelling framework and has therefore not been properly assessed in past studies. This term is non-negligible and potentially in line with the impact from other biomass burning by-products.

We can now proceed to compare the impact of BBA over Amazonia with the effect of ${\mathrm{O}}_{\mathrm{3}}$ on the vegetation that is produced from ${\mathrm{O}}_{\mathrm{3}}$ precursors emitted by forest and grassland fires. Although Pacifico et al. (2015) reported the changes in GPP, their results can be directly compared to the changes in NPP derived from our simulations because the effects of BBA in HadGEM2-ES are predominantly affecting the GPP, whereas the impact on plant respiration is of second order over this region of the world under the present-day climate (Fig. S9). Using the same modelling set-up as in the present study, Pacifico et al. (2015) estimated that present-day ${\mathrm{O}}_{\mathrm{3}}$ produced from precursors emitted by forest and grassland fires in the Amazon region reduces the vegetation GPP by approximately $-\mathrm{230}$ TgC yr${}^{-\mathrm{1}}$ over the same region that was analysed in this study. This is about 2 times, but of the opposite sign, the magnitude of the net impact of BBA estimated in this study (i.e. $+\mathrm{80}$ to $+\mathrm{105}$ TgC yr${}^{-\mathrm{1}}$ for the BBAx1 and BBAx2 scenarios), which includes the climate feedbacks. However, it is important to emphasise that the result from Pacifico et al. (2015) is based on an approach of modelling the ${\mathrm{O}}_{\mathrm{3}}$ effects on photosynthesis that includes a “high” and “low” parameterization for each plant functional type to represent species more sensitive and less sensitive to ${\mathrm{O}}_{\mathrm{3}}$ effects. The $-\mathrm{230}$ TgC yr${}^{-\mathrm{1}}$ decrease in GPP reported there is based on the high sensitivity mode to establish the maximum response. It is also worth noting that due to a lack of knowledge and data on the impacts of ${\mathrm{O}}_{\mathrm{3}}$ on tropical vegetation, the ${\mathrm{O}}_{\mathrm{3}}$ damage parameterization in the work by Pacifico et al. (2015) was derived from data from the temperate and boreal regions. As discussed in the previous paragraph, the BBA-induced DFE is small in our simulations, and if an upper estimate of the BBA were to be considered, it is then possible to argue that BBAs have the potential to virtually counteract the ${\mathrm{O}}_{\mathrm{3}}$ leaf damage resulting from biomass burning in the area. However, while the biomass burning and ${\mathrm{O}}_{\mathrm{3}}$ impacts are potentially of the same magnitude but of the opposite sign, they are not geographically collocated. This means that BBA and ${\mathrm{O}}_{\mathrm{3}}$ do not necessarily affect the same regions of the Amazon rainforest. As reported in Pacifico et al. (2015), ${\mathrm{O}}_{\mathrm{3}}$ tends to show its highest concentrations upwind of the fires which are located over dense areas of broadleaf trees in the model. In contrast to this, the highest AOT from BBA is found downwind of the fires and located over predominantly grassland areas. Future research aimed at assessing the overall net impact of forest and grassland fires on ecosystems through the ${\mathrm{O}}_{\mathrm{3}}$ and DFE effects should therefore consider modelling the two effects simultaneously in a fully coupled framework.

We showed in Sect. 3.3 that the impact of BBA on vegetation over the Amazon rainforest is dominated by the contribution of the term we have referred as climate feedbacks. The (bio)physical mechanisms involved behind this term are numerous, and it is beyond the scope of this paper to completely untangle and quantify them. Future work should seek to understand how aerosol can benefit vegetation productivity when the DFE does not suffice to explain the increase in vegetation NPP. Two working hypotheses for making progress are proposed. First we have noted that BBAs are capable of cooling surface temperatures significantly, which potentially reduces evapotranspiration (ET) and consequently water stress due to a low soil moisture content (Fig. S11a, b). Remarkably, the canopy-level water use efficiency (WUE $=$ GPP$/$ET) is significantly enhanced under higher BBA conditions (Fig. S11d). Given the modest increase in GPP reported earlier, it probably implies that the decline in ET was steeper than the increase in GPP, and this would suggest that vegetation is able to balance water loss and carbon uptake with increasing aerosol concentrations.

Secondly, we suggest that future studies put an emphasis on how BBA can modify the biotic (e.g. rate of carboxylation of the rubisco enzyme, $V_{\mathrm{cmax}}$, leaf temperature) and abiotic factors (air temperature, vapour pressure deficit, PAR, leaf surface temperature, ${\mathrm{CO}}_{\mathrm{2}}$ concentration and air pressure) which control the vegetation response (Lloyd and Farquhar, 2008; Wang et al., 2018). We found that the cooling effect of BBA (Fig. S12a) actually reduces the leaf temperature beyond the $V_{\mathrm{c},\max }$ temperature optimum, which works to reduce plant productivity (Fig. S12c). But the aerosol cooling also lowers the VPD (vapour pressure deficit) which can stimulate stomatal conductance and thus enhance canopy photosynthesis (Fig. S12b). The antagonistic effects from VPD and $V_{\mathrm{cmax}}$ changes are particularly relevant to the sunlit leaves as this population of leaves is mostly rubisco-limited in our modelling framework (not shown). Assessing the role of these ecophysiological mechanisms is critical for developing a better understanding of the ecosystem–climate feedbacks which control the carbon flux from the atmosphere to the land surface, and more attention should be paid to this issue. Further research on the ecosystem–climate feedback will also contribute significantly to understanding the complex relationships between aerosols and ecosystems (e.g. Schiferl and Heald, 2018).

Summary

Intense biomass burning events happen regularly in the vicinity of the Amazon rainforest during the dry season ($\sim$ August–September), releasing huge amounts of trace gases, aerosols, and ozone and aerosol precursors. This potentially leads to very large interactions between chemistry, aerosol, clouds, radiation and the ecosystems.

In this study, we have investigated the impact of biomass burning aerosol (BBA) emissions under present-day conditions on the photosynthesis rate and net primary productivity (NPP) of the Amazon rainforest. Aerosol impacts have many impacts that could influence the ecosystems on a regional scale. Amongst these, light scattering from aerosols is often expected to promote more efficient use of radiation by vegetation through the so-called diffuse PAR fertilisation effect (DFE). To understand the potential impact of BBA in this region, we have implemented an updated representation of plant photosynthesis and carbon uptake that is sensitive to diffuse light radiation in the UK Met Office Hadley Centre HadGEM2-ES Earth system model.

Overall, our simulations indicate that the net impact of BBA increases vegetation NPP by $+\mathrm{80}$ to $+\mathrm{105}$ TgC yr${}^{-\mathrm{1}}$ over the central Amazon Basin (Sect. 3.2). For the first time we have separated the contribution from the individual radiative and climatic processes that contribute to our estimate of the BBA net impact on the vegetation. We found that the increase in diffuse PAR (i) stimulates photosynthesis in the shaded part of the canopy and increases NPP by $+\mathrm{65}$ to $+\mathrm{110}$ TgC yr${}^{-\mathrm{1}}$ in our simulations, (ii) reduces leaf temperature and together with other climatic feedbacks increase NPP by $+\mathrm{67}$ to $+\mathrm{100}$ TgC yr${}^{-\mathrm{1}}$, and (iii) reduces the total amount of radiation, therefore decreasing NPP by $-\mathrm{52}$ to $-\mathrm{105}$ TgC yr${}^{-\mathrm{1}}$, with an overall impact of BBA beneficial for the vegetation.

In our simulations, the DFE from BBA aerosols is small over the analysis region. Our results do not imply, however, that diffuse light is not effective at stimulating vegetation productivity, rather that it is only one of a number of responses to a perturbation in the flux of BBA to the atmosphere. We have discussed some possible reasons why the DFE from BBA appears to be weak in our modelling study (Sects. 3.3 and 4.2). Aerosols are not the only light scatterers present in the atmosphere; clouds too, strongly modify the amount and quality of the radiation reaching the surface. Aerosol-induced DFE impacts may then also depend on cloud cover, which allows for smaller or larger radiative perturbations for the same level of aerosols (e.g. Cohan et al., 2002; Yue and Unger, 2017). Future studies seeking to investigate the DFE of aerosols should therefore critically asses the role played by clouds in providing the baseline diffuse light conditions at the surface before assessing the perturbation associated with aerosol emissions.

The novel result from this study shows that aerosol impacts on vegetation can be significant thanks to the contribution of the climate feedbacks, which are the result from the system adjustment to the aerosol perturbations which ultimately affect vegetation productivity. Those impacts can only be captured when considering the BBA effects in a fully coupled modelling framework. Because the aerosol cooling at the surface has a strong effect on biotic and abiotic processes which control the vegetation response (Wang et al., 2018), future work should invest effort into understanding how the effects of BBA, and other aerosols more generally, can affect the surface energy budget which preconditions photosynthetic activity. This step will certainly become even more relevant as advances in the representation of vegetation physiology and phenology in ESMs are made (e.g. increasing plant functional types or improving vegetation traits), which would likely lead to different vegetation sensitivities to aerosol effects.

Our modelling study specifically aimed at quantifying the changes in the fast ecosystem responses (e.g. NPP or GPP) in response to the effects of BBA. Because the design of our simulations prevents the slow carbon pools to adjust, we cannot investigate how BBA affects carbon allocation, and the potential impact it could have on vegetation structure and dynamics. More research is required to investigate how the impacts of BBA, and indeed all aerosols, on light and on the surface energy budget may alter the onset and shutdown dates of photosynthesis, growing season length, and the canopy structure that provide a feedback to vegetation productivity (Yue et al., 2015). Such feedbacks could become even more relevant under a future warmer climate as anthropogenic aerosol emissions are expected to decrease, while vegetation will continue to experience more and more stressful climatic conditions (e.g. Schiferl and Heald, 2018).

HadGEM2-ES, JULES and SOCRATES codes are available from https://code.metoffice.gov.uk/ (last access: 25 January 2019) for registered users. To register for an account, users should contact their local institutional sponsor. If in doubt, please contact [email protected] for advice stating your affiliate institution and your reason for wanting access.

The MODIS cloud and aerosol products (10.5067/MODIS/MYD06_L2.006, last access: 25 January 2019) (Platnick et al., 2015) are available from https://ladsweb.modaps.eosdis.nasa.gov/ (last access: 25 January 2019). The CERES radiation data are from SSF 1-degree Terra edition 2.8, available from https://ceres.larc.nasa.gov/order_data.php (last access: 25 January 2019) (Smith et al., 2019). GPCP version 2.3 combined precipitation datasets are available from https://www.esrl.noaa.gov/psd/data/gridded/data.gpcp.html (last access: 25 January 2019) (Adler et al., 2018). The FLUXCOM data are available from the data portal of the Max Planck Institute for Biogeochemistry https://www.bgc-jena.mpg.de/geodb/projects/Home.php (last access: 25 January 2019) (Jung et al., 2017b). The CRU datasets are available from http://www.cru.uea.ac.uk/data (last access: 25 January 2019) (Jones et al., 2014). MODIS MOD17A2 NPP product was accessed from https://neo.sci.gsfc.nasa.gov/view.php?datasetId=MOD17A2_M_PSN (last access: 25 January 2019) (Running et al., 2015). The EMDI data are accessible from http://gaim.unh.edu/Structure/Intercomparison/EMDI/ (last access: 25 January 2019) (GAIM, 2019). Figures were prepared using the NCAR Command Language (NCL, version 6.4.0) 10.5065/D6WD3XH5 (UCAR/NCAR/CISL/TDD, 2018).

The supplement related to this article is available online at: https://doi.org/10.5194/acp-19-1301-2019-supplement.

FFM contributed to the text, implemented the parameterisation in the HadGEM2-ES model, set up the HadGEM2-ES, JULES and SOCRATES simulations, and performed the processing and analysis of the model outputs. GAF, LMM, SS and NB contributed to the text, developed the model and implemented the parameterisation in JULES and HadGEM2-ES. PA and JMH contributed to the text, coordinated the SAMBBA project and performed the aircraft deployment.

The authors declare that they have no conflict of interest.

This article is part of the special issue “South AMerican Biomass Burning Analysis (SAMBBA)”. It is not associated with a conference.

Acknowledgements

This work was funded by the Natural Environment Research Council (NERC) South AMerican Biomass Burning Analysis (SAMBBA) project under grant code NE/J010057/1. Gerd A. Folberth was supported by the joint UK BEIS–Defra Met Office Hadley Centre Climate Programme (GA01101) and the European Union's Horizon 2020 Programme for Research and Innovation under grant agreement no. 641816 (CRESCENDO). Florent F. Malavelle and Jim M. Haywood were partly funded by the NERC SWAAMI grant NE/L013886/1. Lina M. Mercado and Nicolas Bellouin were partly supported by the UK Natural Environment Research Council through The UK Earth system modelling project (UKESM, grant no. NE/N017951/1). Paulo Artaxo acknowledges FAPESP (Fundacao de Amaparo a Pesquisa do Estado de São Paulo) projects 2017-17047-0, 2013/05014-0 and 2012/14437-9 and acknowledges the support from LBA Program that is managed by INPA (The Brazilian National Institute for Amazonian Research). Edited by: Hugh Coe Reviewed by: two anonymous referees