2.1 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^∘ by 1.25^∘) with 38 vertical levels with the model top at ∼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_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_L0), 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).

2.2 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^∘ S, 70–53^∘ 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 −100 %, −50 %, 0 %, +100 % and +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^∘ S, 85^∘ W; 15^∘ N, 30^∘ 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. Δ 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).

2.3 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.

Table 1List of model simulations done for the five experiments.

Download Print Version | Download XLSX

2.3.1 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=0.1 µm, and a geometric standard deviation (σ), σ=1.30, are assumed. For aged smoke, r=0.12 µm and σ=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. ∼0.80) absorb more solar radiation than BBAs with higher SSA (e.g. ∼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 ±10 % by scaling the specific scattering (K_sca in m² kg⁻¹) and absorption (K_abs in m² kg⁻¹) coefficients (K_sca in m² kg⁻¹) 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.

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

Download Print Version | Download XLSX

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_{DIFF_OP}, BBAx1_{DIFF_OP} and BBAx2_{DIFF_OP} for the diffuse case and BBAx0_{ABS_OP}, BBAx1_{ABS_OP} and BBAx2_{ABS_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_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.

2.4 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_d, the total PAR, TotPAR, and the climate feedback, clim, such as NPP(f_d, TotPAR, clim). Neglecting the interdependency between the three terms enables the following decomposition:

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_d and TotPAR using the four PAR fluxes previously introduced; one that considers the effect of aerosols (f_d.aer and TotPAR.aer) and one that considers clean-sky conditions (f_d.clean and TotPAR.clean).

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

Download Print Version | Download XLSX

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:

with

where overbars denote quantities averaged over a time period long enough for vegetation fast responses to adjust to the aerosol effects.

Figure 1Global 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.

Download

2.5 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.

Figure 2Multi-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).

Download

3.1 Evaluation

3.1.2 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.

Figure 3Modelled 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 ±1 standard deviation. Dashed lines are the multi-year annual means.

Download

3.1.3 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 −18.0 to −7.5 W m⁻² (i.e. −14.0 to −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 −30 W m⁻², or −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 +6.0 to +12.0 W m⁻² (i.e. approximately +14.0 to +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_d) of +20.0 to +55.0 % (Fig. 3e, f).

Figure 4Showing 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).

Download

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_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_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.

3.2 The net impact of BBA on forest productivity

Figure 4d represents NPP as a function of f_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_d approximately equals 52 %–56 %. The existence of an optimum f_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_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_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_d.

Figure 4c could be used to infer an AOT for which f_d is getting close to the optimum value of 0.55 (Fig. 4d). This would approximately occur at an AOT of ∼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.

Figure 5NPP 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.

Download

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.

Figure 6Mean 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).

Download

Overall, based on the BBAx1 and BBAx2 simulations, we estimate that BBAs increase NPP by about +80 to +105 TgC yr⁻¹, 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⁻¹. 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.

Figure 7Similar 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.

Download

3.3 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_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 +65 to +110 TgC yr⁻¹ 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 −52 to −105 TgC yr⁻¹ 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 +13 and +5 TgC yr⁻¹ 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⁻¹ 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 +67 to +100 TgC yr⁻¹ 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 +41 TgC m⁻¹. The corresponding impact of the reduction in total PAR decreases NPP by −66 TgC m⁻¹. 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 −42 TgC yr⁻¹ 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 (+194 TgC yr⁻¹), resulting in the net impact of BBA on the vegetation to be overall positive (+151 TgC yr⁻¹).

Figure 8Shown on panel (a) the relative changes in NPP (ΔNPP_Net, in grey), the relative changes in NPP due to the change in diffuse fraction (ΔNPP_Frac Diff, in blue), reduction in total PAR (ΔNPP_TOTPAR, in red), the sum of change in diffuse fraction and reduction in total PAR (ΔNPP${}_{f{}_{\mathrm{d}}+\mathrm{TOTPAR}}$, in green, i.e. the DFE), and the climate feedback (ΔNPP_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 (Δ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.

Download

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 ∼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_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.

Figure 9Shown 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_+25 ppm – X0_+25 ppm is the differences between the BBAx1 and BBAx0 simulations with +25 ppm increase in CO₂ concentration).

Download

3.4 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 +28.4 to 38.6 TgC month⁻¹ in August. For the DIFF_OP simulation (respectively ABS_OP) the net change in NPP is +32.1 to 36.2 TgC month⁻¹ (respectively +17.9 to 18.2 TgC month⁻¹). 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⁻¹ 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 CO₂ concentration (experiments 25 and +50 ppm). While increasing atmospheric CO₂ 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 +25 ppm case compared to the main experiment and the +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 CO₂ on the BBA-induced DFE. Note that the atmospheric CO₂ 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 CO₂ concentration should be increased to avoid additional statistical noise produced in the warming climate.