Air quality and health impacts of vegetation and peat fires in Equatorial Asia during 2004–2015

Fire emissions are from FINNpeatSM, described in detail in Kiely et al (2019) and summarised briefly here. FINNpeatSM includes vegetation fire emissions from FINNv1.5 (Wiedinmyer et al 2011). When MODIS fire hotspots are detected on peatland (World Resources Institute 2017) we assume that fires burn into the peat. Emissions from peat fires are estimated from the burn area, peat burn depth, peat density and emission factors (EF). We assume 100 ha of surface burned area for each fire hotspot (as in FINNv1.5), but only 40 ha of peat burn to account for the fact that not all surface fires on peatland will burn into the peat. We estimate the burn depth of the peat based on daily soil moisture from the European Space Agency (ESA CCI SMv04.4) averaged to 2° degree resolution (Liu et al 2012, Dorigo et al 2017, Gruber et al 2017). We assume peat burn depth scales linearly with soil moisture between a maximum burn depth of 37 cm (averaged from Page et al 2002, Usup et al 2004, Ballhorn et al 2009) when soil moisture is low (< 0.15 m³ m⁻³) and a minimum burn depth of 5 cm when soil moisture is high (>0.25 m³ m⁻³). Emission factors (EF) for peat burning are taken as an average of previous studies of burning of Indonesian peat(Christian et al 2003, Hatch et al 2015, Stockwell et al 2016, Nara et al 2017, Wooster et al 2018, Jayarathne et al 2018, Roulston et al 2018). The (EF) for PM_2.5 used for peat fires (22.3 g kg⁻¹) is larger than in other fire emission inventories, such as the Global Fire Emissions Database (GFED4s) and the Global Fire Assimilation System (GFAS) which both use 9.1 g kg⁻¹ (Van Der Werf et al 2010, Kaiser et al 2012).

WRF-chemv3.7.1 was used to simulate PM concentrations across Equatorial Asia (figure 1). The model has been run at 30 km resolution with 33 vertical levels, between the surface and 50 hPa. We used the model to simulate the 6 dry-seasons (August–October) with the greatest fire emissions over the 2004 to 2015 period: 2004, 2006, 2009, 2012, 2014 and 2015. Our domain excludes West Papua, where fires occurred in 2015 (Lohberger et al 2017). All simulations included a 14 day spin up for chemistry at the start of the time period, and with a 24 hour spin up for meteorology every 15–16 day using National Centre Environmental Prediction Global Forecast System (NCEP 2007). In between the meteorology was free running, to allow the model to simulate impacts of fire smoke on meteorology. The MOZART (Model for Ozone and Related Chemical Tracers, version 4; Emmons et al 2010) chemistry scheme was used to calculate gas-phase reactions, with MOSAIC (Model for Simulating Aerosol Interactions and Chemistry; Zaveri et al 2008, Hodzic and Knote 2014) used to represent aerosol processes, separated into 4 bins; 0.039–0.156 µm, 0.156–0.625 µm, 0.625–2.5 µm and 2.5–10 µm. SOA formation from fires in the model is calculated as 4% of the fire emitted CO based on Spracklen et al (2011). A more complete model description can be found in the supplement (table S1 (available online at stacks.iop.org/ERL/15/094054/mmedia)).

Zoom In Zoom Out Reset image size

Anthropogenic emissions are from EDGAR-HTAP2 (Janssens-Maenhout et al 2015) for 2010, biogenic emissions are from MEGAN (Model of Emissions of Gases and Aerosols from Nature; Guenther et al 2006). Following Kiely et al (2019), we inject half of the fire emissions at the surface with the rest spread throughout the boundary layer. For each year, model simulations were completed with and without fire emissions. The contribution of fires to PM concentrations is calculated as the difference between the simulations with and without fire.

Hourly measurements of PM₁₀ (mass concentration of particulate matter < 10 μm aerodynamic diameter) are available from a network of 53 surface sites across Malaysia (Mead et al 2018) for all the periods of this study (figure 1). Hourly PM₁₀ is also available from Pekanbaru in Indonesia for 2013 and 2015, and from Bukit Kototabang in Indonesia for 2004, 2006 and 2009. Weekly averaged PM₁₀ measurements are available from six sites in Indonesia for 2014 and 2015. Hourly measurements of PM_2.5 from 5 locations in Singapore are available for 2014 and 2015, and are averaged to give mean concentrations for Singapore.

Measurements of PM are mainly from urban locations away from the locations of fires. To estimate the PM concentrations from fire at each measurement location we subtract the background PM concentration during months with little fire (months when PM_2.5 fire emissions are <0.1 Tg month⁻¹ across Indonesia).

We averaged hourly data to give daily means, and calculated the fractional bias (FB), Pearson correlation (r), the normalized mean bias factor (NMBF) and normalized mean absolute error factor (NMAEF) (Yu et al 2006) to evaluate the model (supplementary methods).

Population weighted PM_2.5 (PW), a metric of population exposure to PM_2.5 concentrations, was calculated as,

$$\begin{equation*}PW = \,\mathop \sum \nolimits^ \frac{{{C_i}*{P_i}}}{{{P_{tot}}}}\end{equation*}$$

where C_i is the PM_2.5 concentration and P_i is the population of grid cell i, and P_tot is the total population of the domain. The population data is from the Gridded Population of the World, Version 4 (GPWv4) (Center for International Earth Science Information Network and NASA Socioeconomic Data and Applications Center 2016). The total population within the domain is 477 million, with 255 million in the Indonesian part of the domain (total Indonesian population is 263 million).

The long term premature mortality was calculated using the simulated annual mean PM_2.5, with and without fire emissions. PM_2.5 from August from the simulation with no fires was used to represent January to July and November to December. Anthropogenic emissions in the tropics have little seasonal variation, and this method has been used previously to estimate population exposure to fires (Crippa et al 2016, Koplitz et al 2016).

Premature mortality per year, M, from disease j in grid cell i was calculated as,

$$\begin{equation*}{M_{ij}} = {P_i}{I_j}\,\left( {R{R_{jc}} - 1} \right)/R{R_{jc}}\end{equation*}$$

where P_i is the population in i, I_j is the baseline mortality rate (deaths year⁻¹) for j, and RR_jc is the relative risk for j at PM_2.5 concentration, c (µg m⁻³). The baseline mortality rates and the population age composition are from the GBD2017 (Institute for Health Metrics and Evaluation 2019), and the relative risks are taken from the Global Exposure Mortality Model (GEMM) (Burnett et al 2018) for non–accidental mortality (non–communicable disease and lower respiratory infections). The GEMM exposure function was calculated using the relationship between long-term exposure to outdoor PM_2.5 concentrations and mortality, from studies across many countries. The GEMM exposure function was chosen as it incorporates data from a study in China where PM concentrations are regularly high, as is the case in Equatorial Asia. Mean, upper and lower uncertainty intervals from the GEMM have been used to produce mortality estimates with a 95% uncertainty interval. Population count, population age, and baseline mortality rates were kept constant for 2004–2015 to estimate the variation due to changes in exposure only.

To explore differences with previous studies, we also estimate mortality following the method used in Koplitz et al (2016), where the baseline mortality for all causes increases by 1% for every 1 μg m⁻³ increase in annual mean PM_2.5 concentrations below 50 μg m⁻³.

The greatest fire emissions occur between August and October each year, with a secondary peak in January to April (figure 2). The largest dry season emissions occurred in 2015, followed by 2006, 2009 and 2004. All of these years experienced monthly total fire emissions that were greater than 1 standard deviation above the long-term monthly mean. Other years with total dry season emissions above the median were 2012 and 2014.

Zoom In Zoom Out Reset image size

Table 1 compares dry season (August to October) burned area, biomass consumption and emissions for FINNpeatSM and GFED4s inventories (van der Werf et al 2017). Averaged across 2004–2015, FINNpeatSM has a greater burned area compared to GFED4s (fractional bias, FB = 1.01). Dry matter fuel consumption is more comparable (FB = 0.15) due to greater average dry matter consumption per unit area burned in GFED4s (15 189 g m⁻²) compared to FINNpeatSM (6476 g m⁻²), as a result of greater average peat burn depth in GFED4s. Peat makes up half of the average dry matter consumption in GFED, compared to a quarter of the dry matter consumption in FINNpeatSM. The average emissions of CO and CO₂ are similar (FB = −0.04 and FB = 0.07) for the two inventories, while FINNpeatSM has greater dry season PM_2.5 emissions (FB = 0.48) (table 1), due to higher PM_2.5 EF for peat combustion applied in FINNpeatSM (22.3 g kg⁻¹) compared to GFED4s (9.1 g kg⁻¹). The total emissions from fires depends on the percentage of peat burned, as well as the overall dry matter consumption (see supplementary results).

GFED4s uses MODIS burned area (Giglio et al 2013), whereas FINNpeatSM applies a 1 km² burned area to detected hotspots. Previous studies have also found that this method results in FINN having a larger burned area than other emissions inventories in Asia (Vongruang et al 2017), while Liu et al (2020) suggest thick haze in Indonesia in 2015 prevented detection of fires and that MODIS burned area may be underestimated by 93%. In FINNpeatSM, average burn depth is 7.3 ± 3.7 cm, compared to 10.8 ± 4.8 cm in GFED4s. These estimates are lower than many burn depths recorded in the field (Ballhorn et al 2009, Stockwell et al 2016), however field measurements are likely to be taken at large fires where burn depths may be deeper than average (Stockwell et al 2016).

There is a strong correlation between the dry season emissions simulated by FINNpeatSM and GFED4s (r = 0.87–0.98 for different pollutants, supplementary results), although GFED4s emissions have greater interannual variability (figure S1), due to greater variability in peat burn depth (figure S2). Emissions are a product of burned area, burn depth and emissions factors. Compensating differences amongst these variables mean that two emission datasets can predict similar emissions for different reasons. Measurements of burned area, burn depth, and emission factors are needed to help further constrain the emission models.

Figure 3 compares the spatial pattern of average dry season PM_2.5 emissions in FINNpeatSM and GFED4s. In both datasets South Sumatra and Kalimantan are responsible for the majority of fire emissions, with Sumatra accounting for 33%–42% of PM_2.5 emissions and Kalimantan accounting for 52–63%, in agreement with previous studies (Kim et al 2015, Wooster et al 2018).

Zoom In Zoom Out Reset image size

Without fire emissions, the model greatly underestimates PM concentrations across Malaysia and Indonesia (NMBF = −3.72) and the temporal variability across the sites with daily data is poorly simulated (r = 0.27). When fire emissions are included, the model still underestimates observed PM (NMBF = −0.47), although the temporal variability is better simulated (r = 0.51) (figure S3). Most measurements are in urban locations and issues resolving urban-scale pollution are likely to contribute to model underestimation. To overcome this we estimated fire-derived PM from the observations by subtracting measured PM concentrations during periods without fire (see Methods), and compared with the simulated PM concentration from fires (the difference between simulations with and without fires). Figure 4 shows the comparison of simulated and observed fire-derived PM at each site. Across all years, the simulation of fire-derived PM is unbiased (NMBF = 0.14) and the model has reasonable skill in simulating the temporal variability at each site (r = 0.43), although there is year to year and site to site variability (see supplementary results). The NMAEF and FB for the comparison of fires derived PM are also low for each year (NMAEF = 1.07, FB = −0.02; figure S4). Our model skill in comparison against PM₁₀ observations at 52 sites is similar to a previous comparison by Crippa et al (2016) who reported a NMBF of −0.24 for comparison against PM₁₀ observations at two sites in 2015.

Zoom In Zoom Out Reset image size

Table 2 gives the average PM_2.5 concentration across the domain and the population-weighted PM_2.5 exposure for Equatorial Asia due to emissions from fires. PM_2.5 and population-weighted PM_2.5 concentrations are greatest in 2015. In 2004 and 2012 there is greater average population weighted PM_2.5 from fires than for 2009, despite 2004 and 2012 having lower total PM_2.5 fire emissions. This is due to there being more fires in Sumatra in 2012 than in 2009, close to populated areas. Despite having lower emissions than Kalimantan, fires in Sumatra can expose a greater population to poor air quality (Reddington et al 2014, Kim et al 2015, Marlier et al 2015, Koplitz et al 2016). We estimated a population-weighted smoke exposure over July to October of 8.8 µg m⁻³ in 2006 (compared to 8 µg m⁻³ simulated by Koplitz et al (2016)) and 25.6 µg m⁻³ in 2015 (compared to 19 µg m⁻³ by Koplitz et al (2016)).

Fires increase exposure to PM_2.5 concentrations above the WHO recommended limit of 25 µg m⁻³ (World Health Organization 2005) (figure 5). In 2015 fires resulted in an average of 20 million people being exposed to a daily PM_2.5 concentration > 150 µg m⁻³ (figure 5(a)), and 66.5 million people being exposed to daily PM_2.5 concentrations > 25 µg m⁻³ for at least one in two days during August–October (figure 5(b)). Crippa et al (2016) found that 69 million people in Equatorial Asia were exposed to unhealthy air quality for one day in two in 2015, and Mead et al (2018) found that 26 million people in Malaysia were exposed to PM₁₀ levels above the WHO recommended limit of 50 µg m⁻³. For other years we estimate 22.2–51.7 million people were exposed to PM_2.5 concentrations above 25 µg m⁻³ for one day in two (figure 5(b)). The majority of people exposed to poor air quality from fires live in Indonesia (51%–80% of people exposed) and Malaysia (15–30%).

Zoom In Zoom Out Reset image size

Table 2 shows the estimated excess premature mortality, years of life lost, and disability affected life years across the domain resulting from exposure to PM_2.5 from fires. For each year studied, exposure to PM_2.5 from fires resulted in over 13 000 excess premature deaths, 300 000 years of life lost and 500 000 disability affected life years.

The greatest number of excess deaths resulting from fires was in 2015. We estimate exposure to PM_2.5 from fires caused 44 000 excess deaths in 2015, less than 75 600 excess deaths estimated by Crippa et al (2016) or the 100 300 excess deaths estimated by Koplitz et al (2016). This difference is due to different methods of estimating the health impacts of exposure to PM_2.5. Koplitz et al (2016) applied a 1% increase in baseline mortality for all causes of non-accidental death, for every 1 µg m⁻³ increase in annual mean PM_2.5 concentration. When we apply the same function with our simulated PM concentrations we estimate 106 000 premature mortalities in 2015, similar to that estimated by Koplitz et al (2016). In 2006, we estimate exposure to smoke from fires results in 22 100 premature mortalities, greater that the 6 000 excess deaths from cardiovascular mortality estimated by Marlier et al (2012) but less than the 37 600 deaths estimated by Koplitz et al (2016). Using the same relative risk as Koplitz et al (2016), we estimate 42 520 excess premature deaths from the 2006 fires, similar to their estimate. This comparison suggests that the largest uncertainty in health impacts is due to uncertainty in exposure response function (i.e. the sensitivity of health to PM exposure) rather than uncertainty in emissions or PM concentrations. Kushta et al (2018) also found that the majority of uncertainty in long term mortality estimates for Europe is related to the relative risk function. For China, Giani et al (2020) found that the uncertainty in the PM concentrations and the relative risk function contributed similarly to overall uncertainty, increasing the range of estimated mortality when both uncertainties were considered. Similar to this study, they found using a different relative risk function led to a much greater difference in estimated mortality, outside the 95% confidence interval. There may also be mortalities from exposure to fire related air pollution which have not been considered in our study. Jayachandran (2013) suggests that the pollution from the 1997 fires in Indonesia may result in early-life mortality, while we have only calculated health impacts for adults.

Figure 6 shows the regional distribution of excess mortality due to PM_2.5 exposure from fires. The largest mortality occurs in Sumatra, with 38% of the total mortalities due to PM_2.5 exposure from fire. This is due to a large population with close proximity to the fires. Kalimantan, which has a higher proportion of the PM_2.5 emissions than Sumatra (table 1), has an average of 23% of the total mortalities. Averaged across the years, Malaysia accounts for 18% of the mortalities and Singapore accounts for 4%.

Zoom In Zoom Out Reset image size

Figure 7 shows the annual mean population-weighted PM_2.5 and the annual mortality resulting from exposure to PM_2.5 from fires as a function of particulate emission (primary PM_2.5 emissions and SOA formation; see Methods) from fires. For the years we have studied there is a linear relationship between particulate emission and population-weighted PM_2.5 (r = 0.99) and between emission and estimated premature mortality (r = 0.99). For each Tg of particulates emitted from fires, population weighted PM_2.5 increases by 2.1 µg m⁻³, and excess annual premature mortality increases by 2940.

Zoom In Zoom Out Reset image size

A linear relationship between emission and exposure may not be expected; exposure to PM_2.5 and resulting impacts on health depend on the location and magnitude of the emissions, as well as the atmospheric transport of pollution. However, in Equatorial Asia, the location of fires and the direction of pollution transport varies little year to year. Each year, dry season fires occur in similar regions of Equatorial Asia (figure 3), consistent south-easterly winds over South Kalimantan and Sumatra result in similar atmospheric transport patterns (Chang and Wang 2005, Heil et al 2007, Wang et al 2013, Lee et al 2017), and the same areas are exposed to poor air quality (figure S6). This leads to the strong linearity between PM_2.5 emissions, PM exposure, and mortality. The sample size used here is small (n = 6), however, our results indicate that it may be possible to make a simple estimate of PM exposure and health impacts from emissions alone. We used the relationship between PM emission and mortality, to estimate the health impacts from fires across 2004–2015. Total August–October PM_2.5 emissions from 2004 to 2015 were 44.8 Tg, resulting in an estimated 131 700 excess premature mortalities in this period. We note the 6 years studied in detail resulted in a combined total of 113 600 excess premature deaths. We also used this relationship combined with the particulate emission per unit area burned (table 1) to estimate the premature mortality resulting from each 1 km² of land burned. For 2004–2014, we estimate 0.25–0.33 deaths per km² of burned area. For 2015, we estimate 0.58 deaths km⁻², due to the deeper peat burn depth in that year. These numbers provide an indication of the potential magnitude of public health benefits from reductions in fire arising from the moratorium on granting new concession licences for industrial agriculture (Wijedasa et al 2018), peatland restoration (Harrison et al 2019) and fire management (Carmenta et al 2017, Jefferson et al 2020).