Global climate forcing driven by altered BVOC fluxes from 1990 to 2010 land cover change in maritime Southeast Asia

Over the period of 1990–2010, maritime Southeast Asia experienced large-scale land cover changes, including expansion of high-isoprene-emitting oil palm plantations and contraction of low-isoprene-emitting natural forests. The ModelE2-Yale Interactive terrestrial Biosphere global chemistry–climate model is used to quantify the atmospheric composition changes, and for the first time, the associated radiative forcing induced by the land-cover-change-driven biogenic volatile organic compound (BVOC) emission changes (+6.5 TgC y⁻¹ isoprene, −0.5 TgC y⁻¹ monoterpenes). Regionally, surface-level ozone concentrations largely decreased (−3.8 to +0.8 ppbv). The tropical land cover changes occurred in a region of strong convective transport, providing a mechanism for the BVOC perturbations to affect the composition of the upper troposphere. Enhanced concentrations of isoprene and its degradation products are simulated in the upper troposphere, and, on a global-mean basis, land cover change had a stronger impact on ozone in the upper troposphere (+0.5 ppbv) than in the lower troposphere (<0.1 ppbv increase). The positive climate forcing from ozone changes (+9.2 mW m⁻²) was partially offset by a negative forcing (−0.8 mW m⁻²) associated with an enhancement in secondary organic aerosol (SOA). The sign of the net forcing is sensitive to uncertainty in the SOA yield from BVOCs. The global-mean ozone forcing per unit of regional oil palm expansion is +1 mW m⁻² Mha⁻¹. In light of expected continued expansion of oil palm plantations, regional land cover changes may play an increasingly important role in driving future global ozone radiative forcing.

1 Introduction

Recent decades have witnessed large-scale land cover and land-use changes on the maritime continent. More than 4.5 Mha of natural forest were cleared in Indonesia alone over 2000–2010, which is a loss of 4.6 % of year 2000 Indonesian natural forest cover (Margono et al., 2014). Increasing demand for palm oil, produced by oil palm trees (Elaeis guineensis), has simultaneously driven widespread expansion of this agro-industrial tree crop (USDA, 2010). Indonesia and Malaysia cumulatively produce 85 % of the current global palm oil supply (USDA, 2017). The amount of land area planted with oil palm in Indonesia and Malaysia increased from 3.5 Mha in 1990 to 13 Mha by 2010 (Gunarso et al., 2013).

Land cover and land-use changes in Southeast Asia perturb the Earth's system in a variety of ways. Deforestation is a significant threat to Southeast Asian biodiversity (Sodhi et al., 2004), and the land-based carbon emissions associated with forest clearing have greatly contributed to Indonesia's status as one of the world's highest emitters of greenhouse gases (GHGs) (FAO, 2014; WRI, 2015). The magnitude of GHG emissions from deforestation is exacerbated by the pervasiveness of carbon-rich peat soils underlying Southeast Asian tropical forests (Carlson et al., 2012; Hooijer et al., 2010; Page et al., 2002; van der Werf et al., 2009). Peat soil drainage (i.e., drying) promotes oxidation of the sequestered carbon (Miettinen et al., 2017) and is often followed by fire clearing, despite the illegality of this practice in Indonesia (Indrarto et al., 2012). Indonesian forest and peat fires have fueled transnational air pollution episodes (Gaveau et al., 2014; Koplitz et al., 2016), potentially causing more than 100 000 premature mortalities in 2015 (Koplitz et al., 2016).

Above-canopy flux measurements taken in Borneo indicate that, compared to the natural forests of maritime Southeast Asia (MSEA), oil palm plantations are much stronger emitters of the biogenic volatile organic compound (BVOC) isoprene (C₅H₈), with mean midday fluxes about 5 times stronger from oil palm (Langford et al., 2010; Misztal et al., 2011). The simultaneous large-scale contraction of low-isoprene-emitting natural forest area and expansion of high-isoprene-emitting oil palm plantations suggest a land-cover-change-driven increase in regional isoprene emissions over recent decades (Silva et al., 2016; Stavrakou et al., 2014). Measurements indicate that the forests of MSEA emit monoterpenes, a class of BVOCs with chemical formula C₁₀H₁₆, but find negligible monoterpene emissions from oil palm (Langford et al., 2010; Misztal et al., 2011). Both isoprene and monoterpenes are precursors to the short-lived climate pollutant tropospheric ozone (Atkinson and Arey, 2003) and secondary organic aerosols (SOAs) (Carlton et al., 2009; Friedman and Farmer, 2018); as such, perturbations in regional isoprene and monoterpene emissions serve as an additional mechanism by which regional land cover change can affect air quality and climate.

A few studies have used global modeling to investigate the atmospheric composition impacts of Southeast Asian oil palm expansion. Ashworth et al. (2012) analyzed the expected impacts from isoprene emission enhancements associated with the partial replacement of natural forest area with oil palm plantations under a land-use change scenario designed to meet a portion of the projected increase in demand for biofuels in coming years. Warwick et al. (2013) analyzed the impacts from isoprene emission enhancements associated with the total conversion of Borneo from forest to oil palm. Both studies quantified the impacts of the isoprene emission changes first by applying a contemporary NO_x inventory and secondly by assuming increased NO_x emissions near the site of land-use change due to enhanced fertilizer application and increased on-site processing of the palm oil. Based on simulations that apply contemporary NO_x inventories, both studies predict reductions in surface ozone co-located with the isoprene enhancements because the increased VOC serves as a net sink for ozone in the low-NO_x atmosphere (Ashworth et al., 2012; Warwick et al., 2013). When NO_x emission enhancements occur in concert with the land-use-change-driven isoprene emission enhancements, both studies simulate a net increase in ozone production and a concomitant increase in local surface ozone pollution (Ashworth et al., 2012; Warwick et al., 2013). Ashworth et al. (2012) predict local enhancements in annual-mean surface concentrations as high as 11 % for ozone and 10 % for biogenic SOA (as the maximum change in any grid cell) in response to a 5 % increase in regional isoprene emissions and simultaneous enhancements in local NO_x emissions. Warwick et al. (2013) predict local changes in monthly-mean surface ozone as strong +70 % when simultaneously increasing NO_x and isoprene emissions. These studies highlight the potential for significant local surface-level pollution impacts from land-use change.

A recent third study quantified the air quality impacts associated with year 2010 oil palm cover relative to an oil-palm-free scenario (Silva et al., 2016). Using the GEOS-Chem chemical transport model and contemporary inventories for NO_x emissions, Silva et al. (2016) simulate local enhancements in surface pollution as high as 26 % (3–4 ppbv) for ozone and 60 % (about 1 µg m⁻³) for SOA. In Kuala Lumpur, Malaysia, the number of days that register ozone levels higher than the limits recommended by the World Health Organization more than doubled due to regional oil palm expansion (56 days based on 2010 oil palm coverage compared to 23 days in the absence of oil palm; Silva et al., 2016), again highlighting the strong impact of Southeast Asian land cover change on surface pollution.

The year 2000, 16-plant functional type (PFT) land cover distribution map that Silva et al. (2016) apply as the no-oil-palm case is designed for global modeling studies and, therefore, lacks information about the distribution of individual species of vegetation. For example, 7.2 Mha of oil palm existed in 2000 in Indonesia, Malaysia, and Papua New Guinea (Gunarso et al., 2013), yet many global vegetation distributions assign these plantations to one or more of a small number of PFTs. Silva et al. (2016) use a 250 m resolution satellite-based map of year 2010 oil palm plantations to define the contemporary oil palm distribution, and they overlay the oil palm on the palm-free base map, displacing the existing land covers in proportion to their fractional distributions. In reality, the prior land cover of oil palm plantations differs widely by region (Gunarso et al., 2013). Furthermore, the modified land cover distribution from Silva et al. (2016) lacks separate delineations for other pervasive land covers in Southeast Asia, including rubber plantations, which covered at least 6.4 Mha on the maritime continent in 2010 (Gunarso et al., 2013). Rubber trees (Hevea brasiliensis) are very weak emitters of isoprene (Geron et al., 2006; Klinger et al., 2002) but very strong emitters of monoterpenes (Baker et al., 2005; Klinger et al., 2002), which are important SOA precursors (Jokinen et al., 2015). Thus, while the study by Silva et al. (2016) provides evidence of significant surface pollution changes induced by oil palm expansion in Southeast Asia, it provides an incomplete picture of the impact of historical land cover change on atmospheric composition.

This small set of studies focuses on the atmospheric composition impacts induced by altered BVOC emissions from Southeast Asian oil palm expansion, all finding that the downwind impacts are smaller in magnitude than the local impacts near the site of land conversion and isoprene emission changes (Ashworth et al., 2012; Silva et al., 2016; Warwick et al., 2013). Ashworth et al. (2012) speculated a small global forcing impact from the increased isoprene emissions in their land-use change scenario, based on the small simulated global changes in the tropospheric burdens of ozone and the hydroxyl radical (OH). However, no study has directly quantified the global radiative impacts associated with the induced changes in atmospheric composition.

Isoprene perturbations in the tropics may have a particularly powerful impact on longwave radiative forcing (Unger, 2014) because the strong vertical mixing prevalent in the tropics provides a mechanism for surface pollution perturbations to impact ozone concentrations in the upper troposphere (Thompson et al., 1997), where, on a per-molecule basis, ozone changes induce the strongest climate impact due to the thermal contrast with the surface (Lacis et al., 1990). In response to isoprene emission enhancements associated with total conversion of vegetated land to oil palm on the island of Borneo, Warwick et al. (2013) simulate a 20 % increase in ozone at 500 hPa over Borneo and a 20 % increase in peroxyacetyl nitrate (PAN) at 500 hPa downwind of Borneo over the Pacific Ocean. PAN is an organic nitrate that can undergo long-range transport before releasing its reactive NO_x moiety (Moxim et al., 1996), providing a means for ozone formation in remote environments (Kotchenruther et al., 2001). The results of Warwick et al. (2013) suggest that regional isoprene emission changes have the capacity to alter ozone concentrations in the free troposphere and, therefore, induce a radiative forcing.

This study uses the ModelE2-Yale Interactive Terrestrial Biosphere (ModelE2-YIBs) global chemistry–climate model, in conjunction with multiple observational datasets, to quantify the global atmospheric composition changes and, for the first time, the concomitant radiative forcings associated with BVOC emission changes from 1990 to 2010 land cover change in MSEA. The calculations presented here consider changes in emissions of both isoprene and monoterpenes. The applied regional land cover changes are derived from a Landsat-based classification (Gunarso et al., 2013) and account for changes in eight land covers that are prevalent in MSEA, including high-monoterpene-emitting rubber trees and high-isoprene-emitting oil palm trees.

2 Data and methods

2.1 ModelE2-YIBs description

Atmosphere-only simulations employ the NASA GISS ModelE2-YIBs global chemistry–climate model. The YIBs model (Yue and Unger, 2015) is a land surface model embedded in the NASA GISS ModelE2 global chemistry–climate model (Schmidt et al., 2014). The model features 2^∘ latitude × 2.5^∘ longitude horizontal resolution, 40 vertical layers (surface to 0.1 hPa), and a physical and chemical time step of 30 min.

The chemical mechanism includes 156 reactions involving 51 chemical species with full coupling of tropospheric and stratospheric chemistry (Schmidt et al., 2014; Shindell et al., 2006). The troposphere features NO_x-O_x-HO_x-CO-CH₄ chemistry; an explicit representation of isoprene; and a lumped hydrocarbon scheme involving terpenes, peroxyacyl nitrates (PANs), alkyl nitrates, aldehydes, alkenes, and alkanes. The representation of hydrocarbons generally follows Houweling et al. (1998), which is originally derived from the Carbon Bond Mechanism IV (Gery et al., 1989) and the Regional Atmospheric Chemistry Mechanism (RACM; Stockwell et al., 1997), but includes several modifications aimed at representing the wide range of chemical conditions found in Earth's atmosphere, such as the addition of reactions important in low-NO_x conditions, including representation of organic peroxy radical chemistry under low-NO_x conditions and introduction of organic nitrate chemistry. Shindell et al. (2013) describe in detail the recent updates to the tropospheric chemistry scheme, including the incorporation of acetone chemistry (Houweling et al., 1998) and the addition of terpene oxidation (Tsigaridis and Kanakidou, 2007). SOA formation is driven by NO_x-dependent oxidation of emissions of isoprene, monoterpenes, and other reactive VOCs following a volatility-based two-product scheme (Tsigaridis and Kanakidou, 2007). The formation of secondary inorganic aerosols, including sulfate (Bell et al., 2005; Koch et al., 2006) and nitrate (Bauer et al., 2007a), depend on both modeled oxidant levels and the availability of source gases. Primary aerosol types include dust (which provides a surface for heterogeneous chemistry; Bauer and Koch, 2005; Bauer et al., 2007b), black carbon, organic carbon, and sea salt (Koch et al., 2006). Stratospheric chemistry, introduced to the chemical mechanism by Shindell et al. (2006), includes nitrous oxide (N₂O) and halogen (bromine and chlorine) chemistry. Recent updates to stratospheric chemistry are summarized by Shindell et al. (2013) and include changes in the representations of polar stratospheric cloud formation (Hanson and Mauersberger, 1988) and heterogeneous hydrolysis of N₂O₅ on sulfate (Hallquist et al., 2003; Kane et al., 2001). Photolysis rate calculations follow the Fast-J2 scheme of Bian and Prather (2002). At each 30 min time step, the simulated distributions of clouds, ozone, and aerosols are passed to the photolysis code, providing a mechanism for simulated changes in aerosols to impact atmospheric chemistry through modification of photolysis rates (Bian et al., 2003).

Global annual-mean mixing ratios are prescribed for the well-mixed greenhouse gases carbon dioxide (CO₂), methane (CH₄), nitrous oxide (N₂O), and chlorofluorocarbons (CFCs) (Meinshausen et al., 2011; Riahi et al., 2007). Prescribed monthly anthropogenic and biomass burning emissions of reactive gas and primary aerosol species follow the MACCity emissions pathway (Angiola et al., 2010; Granier et al., 2011) for all years, except for 2010, when the interpolated ACCMIP-RCP8.5 dataset (Heil and Schultz, 2014) is applied for biomass burning emissions (as MACCity biomass burning emissions are available only through 2008). MACCity biomass burning emissions were built from the ACCMIP (Lamarque et al., 2010), REanalysis of the TROpospheric chemical composition (RETRO; Schultz et al., 2008), and Global Fire Emissions Database (GFED-v2; van der Werf et al., 2006) datasets (Granier et al., 2011). The interpolated ACCMIP-RCP8.5 emissions were created using simple temporal interpolation of the ACCMIP and RCP8.5 datasets (Heil and Schultz, 2014). Identical surface maps and emission factors were applied in the creation of both the MACCity and interpolated ACCMIP-RCP8.5 biomass burning emissions datasets (Heil and Schultz, 2014). The model dynamically calculates climate-sensitive emissions of reactive compounds for a number of natural sources. These emission sectors, with their primary driving meteorological variables, include: lightning NO_x (moist convection; Price and Rind, 1992; Price et al., 1997), soil NO_x (precipitation and temperature; Yienger and Levy II, 1995), dust (wind speed; Miller et al., 2006), sea salt particles (wind speed; Koch et al., 2006), marine dimethyl sulfide (wind speed; Koch et al., 2006), and the BVOCs isoprene (radiation, temperature, and soil moisture; Arneth et al., 2007; Unger et al., 2013) and monoterpenes (temperature; Lathière et al., 2006), which are described in more detail below.

Leaf-level gas exchange in ModelE2 (Collatz et al., 1991) couples the Farquhar–von Caemmerer kinetic model of photosynthetic CO₂ uptake (Farquhar et al., 1980; Farquhar and von Caemmerer, 1982) to the Ball–Berry model of stomatal conductance (Ball et al., 1987). The environmental inputs used to drive the vegetation biophysics, isoprene emissions, and monoterpene emissions are the values simulated by the general circulation model. The standard YIBs vegetation comprises eight plant functional types: C3-grassland, C4-grassland, crop, deciduous broadleaf forest, evergreen broadleaf forest, evergreen needleleaf forest, shrubland, and tundra. The model code was modified to include four additional land cover types: (1) oil palm plantations, (2) rubber tree plantations, (3) other tree plantations, and (4) dipterocarp evergreen broadleaf forest. Relative to the forests of the Amazon, the tropical forests of MSEA contain a larger proportion of evergreen dipterocarp forests (Hewitt et al., 2010). The dipterocarp evergreen broadleaf forest PFT was added to account for the comparatively lower isoprene emission capacity (Langford et al., 2010; Stavrakou et al., 2014) of the natural forests of MSEA. In each model grid cell, an individual canopy is simulated for each PFT. The canopy radiative transfer scheme divides each canopy into a flexible number of vertical layers (generally 2–16). In each layer, sunlit leaves use direct photosynthetically active radiation (PAR) for leaf-level photosynthesis, while shaded leaves use diffuse PAR (Spitters et al., 1986).

YIBs features a process-based biochemical model of isoprene emission in which the rate of isoprene production dynamically depends on the electron-transport-limited rate of photosynthetic carbon assimilation (Arneth et al., 2007; Unger et al., 2013). For each PFT, the isoprene emission rate depends linearly on the fraction of electrons available to undergo isoprene synthesis, the calculation of which requires prescription of a PFT-specific leaf-level isoprene basal emission rate (BER) at standard conditions. The leaf-level isoprene emission rate additionally depends on the atmospheric CO₂ concentration and the simulated canopy temperature and PAR. The ability of the model to simulate isoprene emissions has been evaluated extensively against multiple above-canopy flux datasets from tropical and temperate ecosystems (Unger et al., 2013). The model simulated the local flux magnitude within a factor of 2 at nine specific measurement sites, some of which correspond to short (weeks-long) measurement campaigns (Unger et al., 2013). Temperature-dependent leaf-level monoterpene emissions, functionally α-pinene, likewise vary by ecosystem type, similarly through prescription of PFT-specific basal emission rates (Guenther et al., 1995; Lathière et al., 2006). Recent work suggests that tropical monoterpene emissions exhibit both a light and temperature dependency (Guenther et al. 2012; Jardine et al., 2015, 2017) that is not included in the emission algorithm here but may be explored in future work.

Aerosol and gas-phase chemistry are fully coupled, and the chemical mechanism is fully coupled to the climate modules (e.g., radiation and dynamics). All simulations apply ozone and aerosol climatologies to the radiation code (Schmidt et al., 2014), but simulated changes in ozone and aerosols can impact calculated photolysis rates. In this study, aerosols do not affect cloud properties. Observation-based, monthly-varying, 5-year-average sea surface temperature and sea ice fields are prescribed according to the Hadley Centre Sea Ice and Sea Surface Temperature dataset (Rayner et al., 2003).

ModelE2 has previously undergone extensive, rigorous validation of simulated present-day tropospheric and stratospheric chemical composition, circulation, and ozone forcing using multiple observational datasets (Shindell et al., 2006, 2013). Shindell et al. (2013) compared simulated monthly zonal-mean total column ozone to that from observations (2000–2010 means) from the Total Ozone Mapping Spectrometer and the Ozone Monitoring Instrument (McPeters et al., 2008), finding that the simulated zonal-mean total column ozone in the tropics shows little bias (<5 %) against measurements for each month, and, in the Northern Hemisphere middle and high latitudes, biases are smaller in the summer months (<10 %) than in the winter months (around 15 %–20 %). Shindell et al. (2013) find only a small negative bias (−0.016 W m⁻²) in the present-day global-average radiative impact of modeled tropospheric ozone relative to TES-derived tropospheric ozone. They note that the strongest biases in ozone concentrations in ModelE2 are generally located in regions where ozone exhibits little effect on radiation (Shindell et al., 2013). More recently, Harper et al. (2018) compared annual-mean ozone concentrations simulated by ModelE2 (representative of year 2005) with an ozonesonde climatology based on measurements taken over 1995–2011 (Tilmes et al., 2012), finding lower model biases at higher pressures (e.g., +2.6 % at 200 hPa compared to +16.9 % at 800 hPa).

A number of uncertainties exist related to both isoprene oxidation chemistry and SOA formation. For example, field measurements, including one campaign over Borneo (Stone et al., 2011), indicate that OH concentrations over the pristine tropical rainforest are much higher than predicted by known isoprene chemistry (Lelieveld et al., 2008; Martinez et al., 2010; Stone et al., 2011). A number of OH-recycling mechanisms associated with isoprene oxidation have been proposed (e.g., da Silva et al., 2010; Lelieveld et al., 2008; Paulot et al., 2009; Peeters et al., 2009). Other researchers have argued that the high measured OH concentrations may be an artifact of the type of instrumentation employed, which results in an artificial inflation of the measured concentrations (Mao et al., 2012). No OH recycling is applied in the OH-initiated isoprene oxidation pathway in the chemical mechanism of ModelE2-YIBs because of the significant uncertainties associated with this aspect of isoprene chemistry. Importantly, Warwick et al. (2013) found that inclusion of OH recycling had little impact on the magnitude or distribution of surface ozone changes induced by biogenic emission changes from the total conversion of Bornean tropical forest to oil palm plantations. However, some studies suggest that the simulation of ozone and other oxidants is sensitive to the isoprene chemical mechanism that is applied (e.g., Archibald et al., 2010; Emmerson and Evans, 2009; Knote et al., 2015). Furthermore, the appropriateness of using α-pinene to represent monoterpenes as a single lumped species in global modeling is an active area of research. Friedman and Farmer (2018) find order-of-magnitude differences in SOA yields for OH oxidation of different monoterpene species, but a clear explanation based on isomer structure remains largely elusive. While Friedman and Farmer (2018) find that the magnitude of the SOA yield from α-pinene is in the “mid-range” of the yields among the analyzed monoterpene species, other studies have shown that this may not be the case for other oxidation pathways (e.g., Draper et al., 2015). The appropriateness of using the two-product scheme for SOA production in global models is likewise an open question (e.g., Tsigaridis et al., 2014) and is discussed in more detail in Sect. 4. Future work would benefit from an exploration of the impact on radiative forcing induced through application of different mechanisms of (1) isoprene photooxidation and (2) SOA formation (e.g., Surratt et al., 2010; Zhang et al., 2018).

2.2 Vegetation datasets

The land cover distributions for MSEA that are applied to the simulations are built from the results of a visual classification of Landsat images that additionally used high-resolution Google Earth images and other datasets (e.g., maps of roads) to aid image interpretation (Gunarso et al., 2013). The Gunarso et al. (2013) classification encompasses 22 land cover types in 30 m × 30 m pixels for the principal oil-palm-producing regions of maritime Southeast Asia – Papua New Guinea, Malaysia (Peninsular Malaysia, Sabah, and Sarawak), and three regions of Indonesia (Kalimantan, Papua, and Sumatra) – nominally for 1990, 2000, 2005, and 2010. The year 2010 Gunarso et al. (2013) dataset relies on some 2009 Landsat images due to intense cloudiness in 2010. As Gunarso et al. (2013) did not perform the classification for Papua New Guinea for 2005, the areal coverage for the grid cells in this region for this year is estimated here through linear interpolation of the year 2000 and 2010 Gunarso et al. (2013) datasets. The Gunarso et al. (2013) classification for 1990 in Indonesia is likewise incomplete; in this dataset, the pixels classified for 1990 are principally those that are oil palm in 1990 or eventually become oil palm. The Gunarso et al. (2013) classification for 2000 was applied to all Indonesian pixels that were not classified in 1990. As such, the 1990 Indonesian land cover map that is applied to the ModelE2-YIBs simulations is a hybrid of 1990 and 2000 data. While Indonesian oil palm cover in 1990 is accurate within the limits of the classification methodology, Indonesian forest cover is presumably lower in the hybrid dataset than it was in reality in 1990, which means that Indonesian forest loss over the period of 1990–2010 is likely underestimated in this study.

Using the 30 m × 30 m pixels from the Gunarso et al. (2013) analysis, including the estimates described above for 2005 in Papua New Guinea and 1990 in Indonesia, the areal coverage of each of the 22 land cover types was quantified for each 2^∘ latitude × 2.5^∘ longitude ModelE2-YIBs grid cell and distributed on the ModelE2-YIBs land surface map. Seventeen of the twenty-two land cover types of the Gunarso et al. (2013) dataset were aggregated into a set of eight cover types (Table S1 in the Supplement). Four of the cover types in the aggregated set (shrubland, crops, C4-grassland, and bare ground) belong to the standard set of YIBs cover types, while the other four cover types (dipterocarp forest, oil palm plantations, rubber plantations, and other tree plantations) are new cover types added to YIBs for this study. The aggregation and subsequent calculation of fractional areal coverage of land for these eight cover types ignored five minor land cover types from the Gunarso et al. (2013) classification: water bodies, coastal fish ponds, mining areas, settlements, and pixels that are unclassified due to clouds.

In the land cover dataset applied to YIBs, 1990–2010 forest loss is likely underestimated because (1) the 1990 Indonesian forest extent is partially derived from the Gunarso et al. (2013) year 2000 map due to unclassified land area in the Gunarso et al. (2013) year 1990 map and (2) deforestation in this dataset represents complete removal of a forest patch and does not account for forest degradation. The ModelE2-YIBs dipterocarp forest PFT combines the “disturbed” and “undisturbed” forest classes from the Gunarso et al. (2013) analysis (Table S1); thus, reduction in forest basal area from the partial logging of a forest does not register as forest loss in the YIBs dataset as long as the affected patch still meets the Gunarso et al. (2013) classification criteria for the disturbed forest class. An analysis of the leaf area index (LAI) of rainforest plots in Central Sulawesi, Indonesia, under different land use regimes found that disturbance of the forest by selective logging reduced the LAI below the 6.2 m² (leaf) m⁻² (ground) value measured for the undisturbed natural forest (Dietz et al., 2007). Removal of “small-diameter” trees reduced LAI to 5.3 m² m⁻², while removal of “large-diameter” trees reduced LAI to 5.0 m² m⁻² (Dietz et al., 2007). Thus, disturbed (selectively logged) forests maintain a high LAI, suggesting that combining the disturbed and undisturbed forest classes into a single PFT is well justified for the purposes of this study.

The Landsat-based YIBs-compatible land cover distributions are applied to the 57 model grid cells covering Peninsular Malaysia, Sumatra, Borneo, and New Guinea (Fig. S1). The simulations apply nonzero areal extents of the four new land cover types only in MSEA. Table 1 shows, for the new land cover types, the assigned physical parameters (including LAI and vegetation height), photosynthetic parameters, and leaf-level basal emission rates of isoprene and monoterpenes.

Table 1 Physical and photosynthetic parameters and isoprene and monoterpene basal emission rates assigned to four new land cover types in the ModelE2-YIBs source code. Plts is plantations.

Parameters: (1) V_cmax25: maximum photosynthetic capacity at 25 ^∘C; (2) I_S: leaf-level isoprene basal emission rate (BER) at standard conditions of incident photosynthetically active radiation (PAR; 1000 µ mol photons m⁻² s⁻¹) and leaf temperature (30 ^∘C); (3) M_S: leaf-level monoterpene BER at standard conditions of incident PAR (1000 µmol photons m⁻² s⁻¹) and leaf temperature (30 ^∘C); (4) SLA: specific leaf area; (5) PAR absorptance: the fraction of PAR photons incident on the leaf that are absorbed by the leaf; (6) vegetation height; and (7) leaf area index (LAI). ^a V_cmax25: value assigned to the standard evergreen broadleaf forest PFT in YIBs. This value is supported by the average of measurements from five common tree species in the forest of Sulawesi, Indonesia (38.4 µmol m⁻² s⁻¹; Rakkibu, 2008). I_S: Upper limit of the BER (<2 µgC g⁻¹ h⁻¹) reported for the species Dipterocarpus obtusifolius (Geron et al., 2006). Three additional tree species in the Dipterocarpaceae family are likewise reported to have low leaf-level isoprene emission rates, where low indicates an emission rate on the order of 0.1 µgC g⁻¹ h⁻¹ (actual numerical rates are not provided; Klinger et al., 2002). M_S: calculated as the mean of the measured leaf-level emission rates for 11 dipterocarp species (Llusia et al., 2014). Height: measurement from a natural forest plot in Malaysian Borneo (Fowler et al., 2011), within the range reported by Dietz et al. (2007) for natural forest plots. LAI: Measurement from natural forest plot in Malaysian Borneo (Fowler et al., 2011). Close to the LAI of 6.2 m² m⁻² estimated for undisturbed forest stands by Dietz et al. (2007). SLA and PAR absorptance: values assigned to the standard evergreen broadleaf forest PFT in YIBs. ^b V_cmax25: corresponds to mature (12-year-old) plantations (Meijide et al., 2017). I_S: Cronn and Nutmagul, 1982; Kesselmeier and Staudt, 1999. M_S: measured emissions for six different monoterpenes (Geron et al., 2006). SLA: average of two measurements (Fan et al., 2015; Legros et al., 2009). PAR absorptance: for measurements of mature leaves (Ritchie and Runcie, 2014). Height and LAI: measurements from 12-year-old commercial plantation in Malaysian Borneo (Misztal et al., 2011). ^c V_cmax25: Kositsup et al., 2009. I_S: Geron et al., 2006; Klinger et al., 2002. M_S: Baker et al., 2005; Klinger et al., 2002. SLA and PAR absorptance: values assigned to the standard evergreen broadleaf forest PFT in YIBs. Height: mean height of rubber trees measured in 49 stands in Peninsular Malaysia (Suratman et al., 2004). LAI: rubber trees are evergreen trees in the humid tropics (Li et al., 2016). Assigned LAI measured for a mature oil palm plantation in Malaysian Borneo (Misztal et al., 2011). ^d Description: this PFT is a combination of the timber plantation and mixed tree crop and agroforest cover types from the Gunarso et al. (2013) land cover classification scheme that is used to build the maritime Southeast Asian land cover distribution maps that are applied to the ModelE2-YIBs simulations. Typical species grown on timber plantations include Gmelina sp., Paraserianthes falcataria, and Acacia mangium (Gunarso et al., 2013). Southeast Asian agroforest plots can contain a diverse array of vegetation, including oil palm, rubber trees, herbaceous crops, and many other tree species used as cash or subsistence crops (e.g., fruit and timber trees) (Scales and Marsden, 2008). Many agroforest plots in Southeast Asia maintain forest-like structural characteristics, despite the fact that they are cultivated, rather than natural, systems (Scales and Marsden, 2008). V_cmax25: value assigned to both the evergreen broadleaf forest and crop classes in YIBs. I_S: Same as dipterocarp forest PFT. Low-isoprene-emitting rubber trees are prevalent in Indonesian agroforest systems (Scales and Marsden, 2008), warranting a low BER for this agroforest-containing PFT. The average I_S (4.7 µgC g⁻¹ h⁻¹) based on two common timber species – Acacia mangium (Klinger et al., 2002) and Gmelina arborea (Singh et al., 2014) – is similar in magnitude to the I_S for the dipterocarp forest PFT, providing further justification for assignment of a low I_S. M_S: same as dipterocarp forest PFT. Similar M_S is reported for Acacia mangium (0.66 µgC g⁻¹ h⁻¹; Klinger et al., 2002). SLA and PAR absorptance: values assigned to the standard evergreen broadleaf forest PFT in YIBs. Height: reported height measured for plots of an agroforest (Dietz et al., 2007) and a timber plantation (Krisnawati et al., 2011). LAI: mean LAI reported for three agroforest plots in Central Sulawesi, Indonesia (Dietz et al., 2007). The timber plantation class from Gunarso et al. (2013) has a partially open canopy, indicating a lower LAI than for the natural forest PFT.

The static land cover distribution applied to the rest of the world is the year 2000 land cover distribution developed for the Community Land Model (CLM; Oleson et al., 2010) using multiple satellite datasets, including retrievals from both MODIS (Hansen et al., 2003) and AVHRR (DeFries et al., 2000). The 16-PFT data were aggregated into the standard set of eight YIBs PFTs (Yue and Unger, 2015). Gridded PFT-specific LAI and vegetation height parameters are prescribed. For the rest-of-world vegetation, prescribed LAIs are derived from the CLM (Oleson et al., 2010); and prescribed heights are the output of a 140-year ModelE2-YIBs simulation (Yue and Unger, 2015) that simulated dynamic carbon allocation, used the same CLM land cover distribution described here, and was forced with year 2000 meteorology from the WFDEI (WATCH Forcing Data methodology applied to ERA-Interim data; Weedon et al., 2014) dataset.

Table 2 Areal extents (Mha) of eight YIBs land cover types in maritime Southeast Asia. Extents encompass only the 57 grid cells that apply the land cover distributions derived from the Gunarso et al. (2013) analysis. The changes in areal extent (in Mha) relative to 1990 are listed in parentheses for 2005 and 2010.

Table 2 shows the areal extents of YIBs land covers in MSEA for 1990, 2005, and 2010. Figure S1 in the Supplement shows the regional land cover distribution for 1990. Over 1990–2010, 11.3 Mha of natural rainforest was lost (−8 % relative to 1990). Forest loss was widespread on Borneo, Sumatra, and Peninsular Malaysia (Fig. S2). Contraction of rubber plantations (−1.4 Mha) was primarily confined to Sumatra and Peninsular Malaysia. The high-isoprene-emitting oil palm class experienced the largest absolute increase in areal extent over the study era (+9.6 Mha, +267 %). Widespread expansion occurred on Sumatra, Borneo, and Peninsular Malaysia.

2.3 Simulation configurations

Table 3 summarizes the configurations of nine global chemistry–climate simulations. Two principal time-slice simulations – 2010land_base and 1990land_base – are used to diagnose the global-mean radiative perturbation associated with the atmospheric composition changes induced by 1990–2010 land cover change in MSEA. The two simulations differ only in terms of the year of the applied maritime Southeast Asian land cover distribution, which is indicated in the simulation names. The regional land cover changes are imposed on a background climate and atmosphere representative of year 2010; that is, both the 2010 land cover distribution in the perturbation simulation and the 1990 land cover distribution in the control simulation are exposed to the climate state, atmospheric CO₂ concentration, and background atmosphere representative of year 2010. The present-day rest-of-world land cover distribution is identical for all simulations. All simulations were run for 13 years, and averages over the final 10 years of output were used for analysis. The impact of 1990–2010 maritime Southeast Asian land cover change, denoted ΔLC, on atmospheric composition and radiative balance is diagnosed as 2010land_base minus 1990land_base (Table 4). Seven additional simulations (Tables 3 and 4) probe the sensitivity of the radiative forcing results to (1) the applied background atmosphere; (2) the degree of regional land cover change; (3) the magnitude of the isoprene BER assigned to the oil palm plantation PFT; and (4) the magnitude of the isoprene BER assigned to the dipterocarp forest PFT.

Table 3 Simulation configurations.

^a The non-land-cover variable boundary conditions include: well-mixed greenhouse gas (CH₄, CO₂, N₂O, and halocarbon) mixing ratios applied to model radiation; CH₄, N₂O, and halocarbon mixing ratios applied to atmospheric chemistry; CO₂ mixing ratio applied to the land biosphere; reactive emissions of carbon monoxide (CO), ammonia (NH₃), NO_x, sulfur dioxide (SO₂), non-methane volatile organic compounds (NMVOCs), black carbon (BC), and organic carbon (OC) from anthropogenic and biomass burning sectors; sea surface temperatures; and sea ice fields. Monthly sea ice and sea surface temperature fields are 5-year averages centered on the simulation year.

Isoprene production in ModelE2-YIBs is calculated as a semi-mechanistic function of photosynthetic carbon assimilation (Unger et al., 2013). Isoprene emissions are sensitive to simulated changes in the parameters that affect photosynthesis, including the background climate state (e.g., temperature, PAR, and soil moisture) and the atmospheric CO₂ concentration. Simulated monoterpene emissions are likewise sensitive to temperature shifts (Lathière et al., 2006). Following emission, the atmospheric processing of isoprene and monoterpenes is influenced by the background atmospheric composition; for example, the availability of NO_x affects the VOC-driven production of both ozone (Sillman, 1999) and SOA (Tsigaridis and Kanakidou, 2007). The calculation 2010land_1990atm minus 1990land_1990atm, denoted ΔLC-1990atm, is designed to test the influence of the prescribed background state upon which the maritime Southeast Asian land cover changes are imposed (Table 4). Relative to the main set of simulations probing Southeast Asian land cover change (i.e., ΔLC), the ΔLC-1990atm simulations prescribe identical changes in Southeast Asian land cover, but impose the changes on a background state representative of year 1990 rather than year 2010 (Table 3). The different background states can lead to different isoprene- and monoterpene-driven impacts on ozone and SOA concentrations from the identical prescribed land cover changes by influencing (1) the magnitude and distribution of the isoprene and monoterpene emission changes associated with the land cover changes and (2) the atmospheric processing of the emitted isoprene and monoterpenes.

A third pair of simulations (2005land_base minus 1990land_base, denoted ΔLC-2005) quantifies the global impacts of 1990–2005 maritime Southeast Asian land cover change (Table 4). Because the land cover distributions in both the control and perturbation simulations experience 2010 background conditions, ΔLC-2005 represents the impacts in 2010 that would have been expected had land cover remained static after 2005. Because land cover in MSEA is rapidly changing and there is uncertainty associated with the classification procedures used to create the land cover maps (Sect. 2.2), such information provides insight into the degree to which the land-cover-change-driven forcing depends on the exact land cover distribution that is applied.

The fourth pair of simulations (2010land_OPber minus 1990land_OPber, denoted ΔLC-OPber) probes the sensitivity of the land-cover-change-induced impacts on atmospheric composition and climate to the magnitude of the isoprene BER applied to oil palm plantations (Table 4). Halving the isoprene BER applied to oil palm plantations brings the simulated 24 h mean isoprene emission rate from oil palm plantations for 2010 close to the rate measured using above-canopy fluxes (described in Sect. 3.3). Similarly, the fifth pair of simulations (2010land_DPTber minus 1990land_DPTber, denoted ΔLC-DPTber) probes the sensitivity of the radiative forcing from 1990 to 2010 land cover change to the magnitude of the isoprene BER applied to dipterocarp forests. In this sensitivity analysis, the dipterocarp forest isoprene BER is increased by a factor of 12, making it equivalent to the isoprene BER assigned to the standard evergreen broadleaf forest PFT in YIBs (Table 1).

3 Results and discussion

3.1 Gross primary production and isoprene and monoterpene emissions

Simulated global gross primary production (GPP) for 2010 is 124 PgC y⁻¹ (simulation 2010land_base), which almost precisely matches an estimate derived from flux-tower measurements that is representative of 1998–2005: 123±8 PgC y⁻¹ (mean ± 1 standard deviation; Beer et al., 2010). The simulated 1990 global GPP of 108 PgC y⁻¹ (simulation 1990land_1990atm) is outside of the 1 standard deviation range of the observation-based mean, but falls within the 95 % confidence interval (102–135 PgC y⁻¹; Beer et al., 2010). Guenther et al. (2012) collated contemporary global annual BVOC emission estimates from the literature, finding ranges of 309–706 TgC y⁻¹ for isoprene and 26–156 TgC y⁻¹ for monoterpenes. The model estimates for 1990 (325 TgC y⁻¹ isoprene and 90 TgC y⁻¹ monoterpenes for simulation 1990land_1990atm) and 2010 (363 TgC y⁻¹ isoprene and 77 TgC y⁻¹ monoterpenes) fall within these ranges. Using the same process-based, leaf-level isoprene production algorithm employed in ModelE2-YIBs, although driven with different forcing datasets, Hantson et al. (2017) predict contemporary isoprene emissions (385 TgC y⁻¹; 1971–2000 mean) that are 18 % higher than those predicted here for 1990 and only 6 % higher than those predicted here for 2010.

Table 4 Calculation methodology. MSEA LCC is maritime Southeast Asian land cover change. The control simulation varies between pairs of simulations.

In MSEA, 2010 isoprene emission rates are generally higher on Sumatra, Peninsular Malaysia, and Borneo than on the island of New Guinea (Fig. S3), which maintains high areal coverage of low-isoprene-emitting dipterocarp rainforests. Rubber plantations contributed 56 % (1.6 TgC y⁻¹) of the regional monoterpene emissions in 2010, while oil palm (8.8 TgC y⁻¹) and shrubs (5.8 TgC y⁻¹) dominated regional isoprene emissions (55 % and 36 %, respectively). The low-isoprene-emitting dipterocarp rainforests, which covered 61 % of the region's land surface in 2010, were responsible for only 8 % of regional isoprene emissions. The strong contributions made by rubber and oil palm plantations to the regional monoterpene and isoprene budgets, respectively, underscore the importance of explicitly accounting for these land covers in regional land use and land cover change analyses.

Regional 1990–2010 land cover change induced a negligible decrease in regional GPP (−0.1 PgC y⁻¹), which is principally attributed to an increase from oil palm (+0.3 PgC y⁻¹) and a decrease from dipterocarp rainforests (−0.4 PgC y⁻¹). Regional land cover change induced annual emissions changes of +6.5 TgC y⁻¹ isoprene and −0.5 TgC y⁻¹ monoterpenes. The land-cover-change-driven net increase in regional isoprene emissions is almost entirely due to expansion of industrial oil palm plantations (+6.4 TgC y⁻¹, +271 % relative to 1990 oil palm isoprene emissions). Regional isoprene emissions from shrubs increased by 4.2 % (+0.2 TgC y⁻¹), associated with a 3.7 % increase in shrubland. The large loss of dipterocarp rainforest had little impact on isoprene emissions (−0.1 TgC y⁻¹), as this PFT is a weak isoprene emitter. Contraction of rubber plantation extent was largely responsible for the reduction in monoterpene emissions (−0.4 TgC y⁻¹).

3.2 Atmospheric composition

Low surface ozone concentrations are simulated for the pristine atmospheres of the tropical forests. In the Southeast Asia study region in 2010 (2010land_base), the less disturbed landscapes of Borneo and New Guinea exhibit lower surface ozone concentrations than the comparatively more industrialized regions of Sumatra and Peninsular Malaysia (Fig. S4). Considering only the grid cells in the 57-grid-cell study area that are majority (>50 %) forest by area, simulated annual-mean surface ozone concentrations in 2010 are 9.0 ppbv in New Guinea (n=10) and 9.5 ppbv in Borneo (n=7). Considering only the grid cells that are >85 % forest by area, simulated ozone concentrations are 7.8 ppbv in New Guinea (n=4) and 7.2 ppbv in Borneo (n=2). Measurements at a rainforest site in Malaysian Borneo in 2008 found daytime surface-level ozone mixing ratios of 5–8 ppbv (Hewitt et al., 2010), providing support for the low ozone concentrations simulated over the highly forested regions of the study area.

1990–2010 land cover change in MSEA drove a reduction in annual-mean surface ozone concentrations over Borneo, Peninsular Malaysia, and Sumatra (Fig. 1). Negligible changes occurred over New Guinea. Isoprene oxidation is more implicated in ozone production and loss rather than SOA formation (whereas monoterpenes are more implicated in SOA formation). Regionally, small surface ozone enhancements are simulated over the marine environment, with maximum enhancements occurring over the ocean to the west of Sumatra. Because this region exhibits low surface ozone concentrations, the small absolute changes (−3.8 to +0.8 ppbv) translate into comparatively large relative changes (−18.3 % to +4.3 %).

Figure 1 Changes in annual-mean ozone mixing ratio (ppbv) due to 1990–2010 maritime Southeast Asian land cover change: (a) Southeast Asian surface-level ozone; (b) global surface-level ozone; (c) tropospheric ozone at 335 hPa; and (d) global-mean ozone profile.

Simulated reductions in surface ozone largely occur in the regions of enhanced isoprene emissions, specifically Peninsular Malaysia, Sumatra, and Borneo (Fig. S5). Surface ozone reductions occurring in response to enhanced VOC emissions in low-NO_x conditions are associated with an increase in the relative importance of ozone destruction reactions (e.g., direct reaction of the VOC with ozone) (Ashworth et al., 2012; Sillman, 1999). NO_x surface emissions in the analysis region were, on average, 0.036 mgN m⁻² h⁻¹ in 2010. Based on atmospheric chemical modeling in conjunction with aircraft and ground measurements in Borneo, NO_x fluxes for 2008 were inferred as 0.009 mgN m⁻² h⁻¹ over a rainforest site and 0.019 mgN m⁻² h⁻¹ over an oil palm plantation (Hewitt et al., 2009). The applied NO_x emissions are slightly higher than, but the same order of magnitude as, the observation-based inferred fluxes, providing support for the “NO_x-sensitive regime” (Sillman, 1999) that is modeled in this study.

Previous studies have shown that the sign and strength of the regional surface ozone response to increased isoprene strongly depend on availability of NO_x (Ashworth et al., 2012; Silva et al., 2016; Warwick et al., 2013). When applying contemporary NO_x emissions from published inventories, both Warwick et al. (2013) and Ashworth et al. (2012) simulated local surface ozone reductions in response to increased isoprene emissions from regional oil palm expansion. The Warwick et al. (2013) study was an idealized simulation that carpeted Borneo in oil palm, while the Ashworth et al. (2012) study applied future projections of oil palm expansion for biofuel production. Both studies found that increasing the NO_x emissions in the region of land conversion (to account for enhanced fertilizer application and industrial processing of the oil palm) enhanced surface ozone concentrations (Ashworth et al., 2012; Warwick et al., 2013). In a study focused on estimating the air quality impacts associated with 2010 oil palm cover compared to a no-oil-palm landscape, Silva et al. (2016) simulated increased surface ozone concentrations over much of the region. They found that some low-NO_x regions (e.g., parts of Borneo) exhibited surface ozone reductions in response to increased isoprene emissions (Silva et al., 2016). The differences in the simulated impacts on regional surface ozone between this study and the Silva et al. (2016) study are likely due to the magnitude of the applied regional NO_x emissions.

The simulated changes in atmospheric composition might be a response not only to altered isoprene and monoterpene emissions, but also to changes in the deposition of atmospheric species induced by changes in leaf density (Wong et al., 2018) or related changes, such as surface roughness, stomatal conductance, and evapotranspiration, that are affected by the applied changes in land cover distribution. Here, the relative changes in regional ozone deposition rates (−19.7 % to +4.3 %) are similar to the relative changes in regional surface-level ozone concentrations (−18.3 % to +4.3 %) from 1990–2010 regional land cover change, in part because the ozone deposition rate depends on the atmospheric concentration change. While increased isoprene emission leading to increased isoprene ozonolysis drives ozone losses near the surface, a formal quantitative attribution analysis disentangling the relative roles of emission and deposition changes requires further complex sensitivity simulations that are beyond the scope of this analysis. In their analysis of Southeast Asian oil palm expansion, Silva et al. (2016) used sensitivity studies to determine that the induced BVOC emission changes, rather than altered deposition rates from LAI changes, were almost exclusively responsible for the simulated surface ozone changes.

While the strongest impacts of regional land cover change on annual-mean surface ozone are confined to Southeast Asia, weak long-range enhancements are simulated, particularly over the tropical ocean (<1 ppbv over Indian Ocean, <0.5 ppbv elsewhere) (Fig. 1). Enhanced isoprene oxidation drives an increase in the formation of alkyl nitrates, which are NO_x reservoirs. Alkyl nitrates sequester reactive NO_x upon formation and can undergo long-range transport, eventually releasing the NO_x far from the source region (Atherton, 1989). Thus, alkyl nitrate perturbations provide a mechanism for ozone perturbations to occur far from the site of the hydrocarbon emission change (Hollaway et al., 2017), here, resulting in ozone enhancements over the tropical ocean.

With decreasing atmospheric pressure, the long-range impact on ozone spreads beyond the tropics and generally grows in magnitude (Fig. 1). The global-mean ozone enhancement increases in magnitude as altitude increases from the surface to ∼100 hPa (Fig. 1). Considering the troposphere, the global-mean ozone enhancement from regional land cover change is on the order of 0.5 ppbv in the upper troposphere (e.g., at 237 hPa), compared to <0.1 ppbv in the lower troposphere (at pressures >875 hPa). The maximum relative change in the global-mean ozone mixing ratio is +0.6 %, occurring at 160 hPa. While enhanced alkyl nitrate formation can likewise play a role in the free tropospheric ozone changes (as for surface ozone), additional mechanisms of ozone change (including in situ production from transported isoprene and its degradation products) are also implicated, as described below.

The land-cover-change-driven maritime Southeast Asian isoprene perturbations occur in a region of deep convection (Folkins et al., 1997). The co-location of emission perturbations and strong vertical transport results in changes in chemical composition throughout the atmosphere (Fig. 2). The strongest zonal-mean enhancements in isoprene occur in the lower troposphere at pressures >850 hPa (maximum $=+\mathrm{88}$ %, occurring near the equator). As a doubly unsaturated hydrocarbon, isoprene is highly reactive in the oxidizing atmosphere. Most of the increased isoprene is rapidly oxidized in the lower atmosphere, yet some is transported to the middle and upper troposphere by strong tropical convection. The isoprene enhancements are generally weaker in the middle and upper troposphere relative to the lower troposphere; considering pressures <850 hPa, the relative change peaks at +37 % at 139 hPa, 5^∘ N.

Upper-tropospheric enhancements in isoprene have been predicted by global model simulations (Collins et al., 1999) and have been observed in convective events over Canada (Apel et al., 2012). When transport-driven isoprene enhancements occur in the presence of lightning NO_x emissions, upper-tropospheric ozone enhancements are predicted (Apel et al., 2012). Lightning NO_x is prevalent in the tropical upper troposphere (Beirle et al., 2010; Boersma et al., 2005), and annual-mean upper-tropospheric ozone enhancements are simulated in this study. The largest zonal-mean tropospheric ozone changes induced by regional land cover change occur in the tropical near-tropopause region (Fig. 2), which is particularly important from a climate perspective because the climate impact of a unit of tropospheric ozone change increases with altitude and is maximized at the tropopause (Lacis et al., 1990). The maximum zonal-mean tropospheric ozone enhancement is +1.4 ppbv (99 hPa, 5^∘ S).

Figure 2 Changes in annual zonal-mean concentrations of (a) isoprene, (b) ozone, (c) formaldehyde, (d) CO, (e) alkyl nitrates, (f) NO_x, and (g) HO_x from 1990–2010 maritime Southeast Asian land cover change. Black traces indicate tropopause location.

Formaldehyde (HCHO) is a high-yield oxidation product of isoprene (Wolfe et al., 2016). In response to maritime Southeast Asian land cover change, enhanced zonal-mean formaldehyde mixing ratios are simulated along the equator from the surface to the upper troposphere (Fig. 2). The strongest changes occur in the lower troposphere, with weaker enhancements in the middle and upper troposphere. Increased upper-tropospheric formaldehyde can result from direct convection of formaldehyde or from in situ production following convection of its precursors (e.g., isoprene).

Laboratory measurements suggest that, in high-NO_x conditions, 60 % of isoprene-derived carbon in the atmosphere eventually becomes CO (Miyoshi et al., 1994). Zonal-mean CO enhancements are simulated for most of the atmosphere, with peak upper-tropospheric enhancements >1 ppbv (Fig. 2). CO is itself an ozone precursor, contributing to the simulated ozone enhancements. The isoprene and formaldehyde change signatures suggest vertical transport in the absence of significant horizontal mixing, which is expected given their short lifetimes on the order of hours. With a lifetime on the order of a few months, CO can undergo significant mixing before being destroyed, explaining why its perturbations are more widely distributed throughout the atmosphere.

In OH-initiated oxidation of isoprene in the presence of NO_x, isoprene nitrates can form from direct reaction of isoprene peroxy radicals with nitric oxide (NO) (Lockwood et al., 2010). In ModelE2-YIBs, 15 % of isoprene + OH reactions in the presence of NO_x produce a primary alkyl nitrate molecule (alkyl nitrates can also be formed from other oxidation pathways) (Shindell et al., 2003). Peak zonal-mean alkyl nitrate enhancements are simulated as +19 % in the lower troposphere, +5 % in the upper troposphere, and +4 % in the stratosphere (Fig. 2). Cold temperatures extend the lifetime of alkyl nitrates in the upper troposphere (Apel et al., 2012), resulting in alkyl nitrate perturbations that are more strongly mixed throughout the upper troposphere and lower stratosphere relative to the lower and middle troposphere (Fig. 2). Sequestration of the reactive NO_x into the alkyl nitrates contributes to the simulated reduction in NO_x mixing ratios in the upper troposphere and lower stratosphere (Fig. 2). Enhanced alkyl nitrate formation is also responsible for surface ozone enhancements over the tropical ocean, as described above.

The land-cover-change-driven isoprene perturbations likewise impact the concentrations of the HO_x (OH+HO₂) radical family. Collins et al. (1999) found that enhanced convection of isoprene and its oxidation products to the upper troposphere enhances HO_x production, principally through increased photolysis of the oxidation products, including formaldehyde. Here, annual zonal-mean HO_x concentrations are enhanced in the tropics from the surface to the upper troposphere in a pattern that largely mimics that of the formaldehyde perturbations (Fig. 2).

Cumulatively, these changes demonstrate that land cover change in MSEA has the capacity to alter the chemical composition of the upper troposphere. Increased isoprene emissions in a region of deep convection results in enhanced upper-tropospheric concentrations of isoprene and its degradation products, which contribute to enhanced ozone mixing ratios in the climatically important near-tropopause region.

3.3 Radiative forcing

Instantaneous radiative forcings are reported as mean ±1 standard deviation, calculated over 10 model years. The uncertainties represent internal model variability. In this study, the quantified radiative perturbations arise specifically from ozone and aerosol changes driven by regional land cover change. Only direct aerosol–radiation interactions are considered for aerosols.

Table 5 Global annual-mean radiative forcing (mW m⁻²) from changes in atmospheric composition induced by maritime Southeast Asian land cover change. Mean ± 1 standard deviation, calculated over 10 model years.

The atmospheric composition changes induced by 1990–2010 maritime Southeast Asian land cover change resulted in a positive annual global-mean radiative forcing of $+\mathrm{8.4}\pm \mathrm{0.7}$ mW m⁻² (Table 5). The global ozone perturbation induced a positive forcing of $+\mathrm{9.2}\pm \mathrm{0.7}$ mW m⁻², offset only slightly by a negative forcing ($-\mathrm{0.8}\pm \mathrm{0.1}$ mW m⁻²) induced by a 1.4 % enhancement (+6.5 Gg) in the global burden of largely reflective SOA particles. (The regional change in SOA is plotted in Fig. S6). The ozone radiative forcing distribution (Fig. 3; range: −10.4 to +37.6 mW m⁻²) largely reflects the pattern of mid- and upper-tropospheric ozone changes (Fig. 1), such that the magnitude of the positive forcing largely tracks the magnitude of the ozone enhancement, particularly in the tropics. On a per-molecule basis, enhancements in ozone have a stronger impact on longwave forcing the nearer they are to the tropopause (Lacis et al., 1990), which explains the close association between the patterns of upper tropospheric ozone changes and the ozone forcing magnitude. The strongest ozone forcings occur over the tropical Indian Ocean.

Figure 3 Annual-mean ozone forcing (mW m⁻²) due to 1990–2010 maritime Southeast Asian land cover change.

The sensitivity studies investigate the uncertainty in the forcing magnitude. For all five drivers, the negative SOA forcing is much smaller in magnitude than the positive ozone forcing, resulting in a total forcing that is positive in sign (Table 5). This consistency provides confidence in the positive sign of the total forcing induced by land-cover-change-driven perturbations in atmospheric composition. Of all of the sensitivity analyses, ΔLC-OPber results in the smallest magnitude forcing for both ozone and SOA (Table 5). The oil palm isoprene BER applied to the simulations used for the ΔLC-OPber calculations was half the magnitude of the published leaf-level BER (Table 1) applied to the simulations used for the ΔLC calculations.

Application of the halved BER brings the simulated 24 h mean isoprene emission rate from oil palm plantations for 2010 (3.8 mgC m⁻² h⁻¹ in simulation 2010land_OPber vs. 7.6 mgC m⁻² h⁻¹ in simulation 2010land_base) in line with the rate measured in 2008 above an industrial oil palm estate in Malaysian Borneo (3.9 mgC m⁻² h⁻¹; Misztal et al., 2011). However, the measured flux, based on a short 12-day measurement window, is sensitive to the prevailing meteorological conditions (Misztal et al., 2011). Misztal et al. (2011) note that the isoprene flux would have been twice as high under the meteorological conditions (PAR and temperature) that they measured during the period prior to initiation of their flux measurements; such a rate would match the rate calculated for the simulation that applied the published leaf-level BER (2010land_base). Thus, the isoprene flux measurements provide confidence in the order of magnitude of the simulated isoprene emissions from oil palm plantations. The ozone forcing calculated using the ΔLC-OPber sensitivity analysis (+4.0 mW m⁻²) is considered to be the lower bound of a best estimate of the global-mean ozone forcing from maritime Southeast Asian land cover change.

Increasing the isoprene BER for the dipterocarp forest PFT by a factor of 12 increases the magnitude of the 1990–2010 isoprene emissions reduction associated with the large contraction in areal cover of this PFT (−1.3 TgC y⁻¹ for ΔLC-DPTber vs. −0.1 TgC y⁻¹ for ΔLC). This decreases the magnitude of the net enhancement in isoprene emissions from total land cover change in the study region (+5.3 TgC y⁻¹ for ΔLC-DPTber vs. +6.5 TgC y⁻¹ for ΔLC). The positive ozone forcing for ΔLC-DPTber is still two-thirds of the magnitude of the ozone forcing for ΔLC (Table 5), despite the factor-of-12 enhancement in the assigned dipterocarp forest BER. Using the 2010 simulation that applies the dipterocarp forest isoprene BER taken from leaf-level measurements (2010land_base simulation), the simulated 24 h mean emission rate of isoprene from maritime Southeast Asian rainforests for 2010 (0.12 mgC m⁻² h⁻¹) is one-third of the eddy-covariance-based measured emission rate for 2008 from a natural forest in Malaysian Borneo (0.35 mgC m⁻² h⁻¹; Langford et al., 2010). The measured flux was shown to be highly sensitive to the meteorological conditions (e.g., wet season flux is 0.22 mgC m⁻² h⁻¹, dry season flux is 0.47 mgC m⁻² h⁻¹; Langford et al., 2010). The results of ΔLC-DPTber suggest that increasing the dipterocarp forest BER by a factor of roughly 3 (to force alignment of the simulated isoprene emission magnitude with the 2008 measurements) would have little impact on the total forcing magnitude (because increasing the BER by a factor of 12 had only a slight impact on the forcing). Because the isoprene emissions capacity of oil palms is so strong relative to that of the dipterocarp forest, the isoprene emissions changes from 1990 to 2010 regional land cover change are dominated by the oil palm PFT. The total land-cover-change-driven forcing is more sensitive to uncertainty in the magnitude of the oil palm BER than to uncertainty in the magnitude of the dipterocarp forest BER. The limited importance of the isoprene emission changes from dipterocarp forest loss additionally indicates that the potential underestimate in forest loss in the land cover dataset applied to these simulations (Sect. 2.2) is unlikely to have a strong impact on the magnitude of the quantified forcing.

The total forcing associated with 1990–2005 land cover change (ΔLC-2005) is 69 % of the forcing associated with 1990–2010 land cover change (ΔLC), indicating that 31 % of the total 1990–2010 forcing is associated with land cover change that occurred over the short 2005–2010 period. The ΔLC-2005 analysis investigates the sensitivity of the simulated forcing to uncertainty in the land cover distribution maps. The ozone forcing is tightly linked to changes in the areal cover of oil palm, and the global-mean ozone forcing per Mha of regional expansion of oil palm plantations is +1 mW m⁻² Mha⁻¹ (the same value is obtained regardless of whether the output from the ΔLC analysis or the ΔLC-2005 analysis is used).

The land cover maps used in this study are based on the land cover classification of Gunarso et al. (2013), which identified large (>1000 ha) patches of oil palm cover. Estimates suggest that around 40 % of Indonesian oil palm area in 2010 (and 26 % in 1990) was associated with smallholders, in contrast to state-owned or private companies (Indonesian Ministry of Agriculture, 2017; Lee et al., 2014). In Malaysia, the smallholder estimate is likewise around 40 % (Vermeulen and Goad, 2006, citing the Malaysian Palm Oil Board). Smallholder plantings are either schemed (contiguous with a larger estate) or independent (Vermeulen and Goad, 2006). In the latter case, plantations are typically on the order of 2–50 ha (Vermeulen and Goad, 2006), far smaller than the patches identified by the Gunarso et al. (2013) analysis. Gunarso et al. (2013) note that it is unknown what proportion of the total smallholder plantation area is excluded from the total oil palm areal cover quantified in their analysis. Nonetheless, it is likely that the oil palm areal extent from the Gunarso et al. (2013) remote sensing analysis is an underestimate of the true areal coverage of oil palm in maritime Southeast Asia.

The 1990–2010 change in oil palm cover dominates the land-cover-change-driven isoprene emissions changes in this region (Sect. 3.1); thus, the underestimate in oil palm areal cover likewise represents an underestimate in the global-mean ozone forcing from regional land cover change. Taking into account the smallholder estimates for Indonesia and Malaysia, the total regional expansion of oil palm cover for 1990–2010 increases to +16 Mha, which is considered to be an upper bound. Applying the land-area-based global-mean ozone forcing for regional oil palm expansion (+1 mW m⁻² Mha⁻¹), based on the ΔLC and ΔLC-2005 analyses, gives an estimate of +16 mW m⁻² for the global-mean ozone forcing from regional land cover change. This value is considered to be the upper bound of a best estimate of the global-mean ozone forcing from maritime Southeast Asian land cover change.

There is little variability in the magnitude of the forcing associated with the prescribed background atmosphere as application of the year 1990 background atmosphere (ΔLC-1990atm) in lieu of the year 2010 background atmosphere (ΔLC) had little impact on the ozone and SOA forcings induced by land cover change (Table 5). The fact that local surface ozone reductions, rather than enhancements, are simulated in response to enhanced isoprene emissions (Fig. 1) suggests a NO_x-sensitive chemical environment.

Taking into account the results of the sensitivity simulations, the best estimate of the global-mean forcing from ozone changes induced by regional 1990–2010 land cover change is +9 mW m⁻², with a range of +4 to +16 mW m⁻². The quantified range accounts only for uncertainties probed by the sensitivity studies.

4 Conclusions and future work

The best estimate of global-mean forcing from isoprene and monoterpene emission perturbations driven by regional land cover change – quantified here using simulations that apply satellite-derived land cover distributions and measured leaf-level isoprene and monoterpene BERs – indicates a positive forcing ($+\mathrm{8.4}\pm \mathrm{0.7}$ mW m⁻²), which is a warming impact. In absolute terms, the quantified forcing from 1990 to 2010 maritime Southeast Asian land cover change is small, particularly in comparison to the forcing associated with industrial-era perturbations of well-mixed greenhouse gases (e.g., Myhre et al., 2013). However, the ozone perturbations associated with changes in global anthropogenic emissions of non-methane VOCs over the industrial era (1750–2011) induced a global-mean stratospherically adjusted forcing on the order of +30 mW m⁻² (Myhre et al., 2013), which is only around 3 times the magnitude of the instantaneous ozone forcing associated with 1990–2010 land cover change in MSEA (+9.2 mW m⁻²). For comparison, the global ozone forcing driven by the 1990–2010 land cover change in MSEA is at the low end of the range of estimates for ozone forcing from global anthropogenic emission source sectors in year 2000 (+5 to +80 mW m⁻²): for example, industry is +15 mW m⁻²; household biofuel is +28 mW m⁻²; road transport is +50 mW m⁻²; power is +53 mW m⁻²; and biomass burning is +71 mW m⁻² (Fuglestvedt et al., 2008; Unger et al., 2010). A multi-model study found that 20 % reductions in NMVOCs (about 2–4 TgC y⁻¹) in four large world regions (North America, East Asia, Europe, and South Asia) in 2001 led to global ozone forcings around −1 mW m⁻² (Fry et al., 2012).

The climate forcing quantified here represents the forcing induced by atmospheric composition changes driven by 1990–2010 land cover change in MSEA. Roughly 27 % of the 2010 oil palm plantation areal cover in this region already existed in 1990 (Table 2). Applying the land-area-based global-mean ozone forcing for regional oil palm expansion that was calculated here (+1 mW m⁻² Mha⁻¹), regional oil palm expansion over the modern era is responsible for a global-mean forcing of +12.7 mW m⁻² from induced ozone changes. Based on government-awarded leases, Carlson et al. (2012) project that at least 9.4 Mha of additional land in Indonesian Borneo alone will be converted to oil palm plantations over 2010–2020, indicating that additional climate forcing is expected from regional land cover change in coming years. For comparison, oil palm expansion over 1990–2010 for the entirety of MSEA was +9.6 Mha (Table 2).

The sensitivity analyses indicate that important factors driving uncertainty in the forcing include (1) uncertainty in the magnitude of the isoprene BER for oil palm and (2) uncertainty in the areal extent of oil palm expansion. The simulations find that the expansion of oil palm plantations is responsible for almost all of the land-cover-change-driven net increase in regional isoprene emissions (Sect. 3.1). Because the total magnitude of isoprene emissions from oil palm plantations scales linearly with the assigned leaf-level BER, the atmospheric composition changes and concomitant forcing from regional land cover change are strongly dependent on the magnitude of the assigned isoprene BER for oil palm plantations. An improved understanding of the isoprene BER for oil palm would strongly benefit the effort to quantify the environmental impacts of the recent large-scale changes in land cover in this sensitive region.

Our study has several limitations. The radiative forcing results are likely sensitive to the isoprene chemical mechanism, SOA production scheme, and convective transport and atmospheric transport schemes in the model. For example, this study applies the two-product scheme for SOA production (Tsigaridis and Kanakidou, 2007), but the appropriateness of using such schemes in global models is still under debate (e.g., Tsigaridis et al., 2014). Many recent global SOA model studies use fixed SOA yields for calculating SOA production from isoprene and monoterpene oxidation (e.g., Rap et al., 2018; Scott et al., 2017, 2018). For the ΔLC analysis, the global-mean SOA radiative forcing per unit of SOA burden change is −115 mW m⁻² Tg⁻¹. This value is largely consistent across the sensitivity analyses, ranging from −112 to −119 mW m⁻² Tg⁻¹. This metric can be used to estimate the SOA radiative forcing induced by the simulated isoprene and monoterpene emission changes under the assumption of fixed SOA yields. Assuming fixed SOA yields of 10 % from the simulated monoterpene emission changes (e.g., Tsigaridis et al., 2014) and 1 % from the simulated isoprene emission changes (lower end of range suggested by Kroll et al., 2005), in conjunction with the SOA forcing per burden metric, results in an SOA forcing of −2.5 mW m⁻² from 1990 to 2010 land cover change (i.e., ΔLC analysis). The SOA radiative forcing based on fixed SOA yields is more than 3 times stronger than, but of the same sign as, the SOA radiative forcing calculated by the global model; in both cases, the SOA radiative forcing is negligible and partially offsets the positive forcing from ozone. For the ΔLC analysis, the cumulative radiative forcing, considering impacts of both ozone and SOA changes, is +8.4 mW m⁻² computed by the model and +6.7 mW m⁻² computed using the simulated ozone forcing plus the SOA forcing computed here using fixed SOA yields. That is, using fixed SOA yields, the total radiative forcing would be slightly smaller in magnitude than, but the same sign as, the forcing simulated by the model. Several recent studies have applied slightly larger SOA yields: +14.3 % from monoterpenes and +3.3 % from isoprene (by mass; Rap et al., 2018; Scott et al., 2017, 2018). Using these larger SOA yields for the ΔLC analysis results in an SOA forcing of −19.4 mW m⁻² and a total radiative forcing, taking into account the ozone forcing, of −10.2 mW m⁻², which is the opposite sign of that simulated by the model (+8.4 mW m⁻²). This analysis indicates that uncertainty associated with biogenic SOA yields from isoprene and monoterpene oxidation has a strong influence on the quantified forcing.

Seasonal variation in isoprene BERs has been observed for some tree species (e.g., Geron et al., 2000). Similarly, based on above-canopy isoprene flux measurements, observed meteorological variables, and the isoprene emission algorithms of the Model of Emissions of Gases and Aerosols from Nature (MEGAN; Guenther et al., 2006) model, Hewitt et al. (2011) suggest that isoprene emissions from oil palm plantations are under circadian control at the canopy scale; that is, the isoprene BER exhibits diurnal variability, peaking in the early afternoon. Keenan and Niinemets (2012) subsequently argued that the apparent circadian control calculated by Hewitt et al. (2011) is likely caused by inappropriate assignment of model parameters in the MEGAN model. One study reports apparent circadian control of oil palm isoprene emissions at the leaf level (Wilkinson et al., 2006). Like other global models, the YIBs model supports application of only a single invariant isoprene BER for each PFT or land cover type and, therefore, might misrepresent isoprene emission magnitudes if oil palm isoprene emissions are, in fact, under circadian control.

Previous work has shown that land cover change in MSEA impacts the environment in numerous ways, including by threatening plant and animal diversity (Sodhi et al., 2004), driving large GHG emissions (FAO, 2014; WRI, 2015), and damaging air quality through vegetation and peat burning (Gaveau et al., 2014; Koplitz et al., 2016). The simulations presented here indicate that isoprene and monoterpene emission perturbations provide an additional mechanism by which regional land cover change impacts the environment. While the impact on global radiative forcing is small, the ozone radiative forcing exceeds +37 mW m⁻² in some localities. This forcing mechanism can be expected to grow in importance in future years if oil palm expansion continues at a high rate.