Climate and health implications of future aerosol emission scenarios

We used a combination of two global climate models to calculate the effects of different aerosol emission scenarios on climate and air quality. To simulate aerosol radiative forcing and aerosol concentrations, we used the aerosol-climate model ECHAM-HAMMOZ (ECHAM6-HAM2.1) (Zhang et al 2012, Stevens et al 2013). ECHAM-HAMMOZ is a general circulation model that solves atmospheric circulation and physics at a horizontal resolution of T63 (roughly 1.9°×1.9°) and accounts for emissions, transport, removal, and microphysics of aerosol particles. The aerosol module includes aerosol species of sulfate, organic carbon (OC), black carbon (BC), sea salt, and mineral dust. Excluding nitrate, ammonium, and secondary organic aerosols means that anthropogenic aerosol forcing and PM_2.5 contribution are likely lower than they would be if these species were included (e.g. Scott et al 2014, Xu and Penner 2012). The model includes both aerosol direct and indirect effects (Lohmann et al 2007) with a semi-empiric cloud droplet activation scheme (Lin and Leaitch 1997). We conducted the ECHAM-HAMMOZ simulations with nudging of the model meteorology towards ERA-Interim reanalysis data (Dee et al 2011) and prescribed monthly climatological (1979–2008) sea surface temperature (Hurrell et al 2008). We used aerosol emissions from the ACCMIP Emission database (http://accmip-emis.iek.fz-juelich.de/data/accmip/gridded_netcdf/accmip_interpolated/T63/) for the historical (Lamarque et al 2010) and RCP (van Vuuren et al 2011a) periods. Natural sea salt, dust, and dimethyl sulfide emissions (DMS) were calculated online. Sea salt emissions were based on a combination of source functions by Monahan et al (1986) and Smith and Harrison (1998) (see Guelle et al (2001) for details), dust emissions were calculated according to Tegen et al (2002) with East Asia soil properties adopted from Cheng et al (2008) (see Zhang et al (2012) for details), and DMS emissions were based on Nightingale et al (2000). For a detailed evaluation of the aerosol concentrations in ECHAM-HAMMOZ, we refer the reader to Zhang et al (2012) and Partanen et al (2013).

ECHAM-HAMMOZ is computationally very expensive and cannot simulate long-term impacts on climate owing to its reliance on prescribed sea surface temperatures. Therefore, we used the University of Victoria Earth System Climate Model (UVic ESCM version 2.9) (Weaver et al 2001, Eby et al 2009) with aerosol forcing from the ECHAM-HAMMOZ runs to assess the wider climate implications of different aerosol emission pathways. The UVic ESCM is an intermediate complexity earth system model that includes a one-layer atmosphere with monthly-varying prescribed winds and simple energy-moisture balance equations, three-dimensional ocean circulation with 19 vertical levels, sea ice, terrestrial vegetation, and ocean biogeochemistry. The model is thus capable of simulating long-term effects on climate and carbon cycle, and its results are comparable to more complex earth system models (Eby et al 2013, Joos et al 2013, Arora et al 2013).

Our main interest in this study was the effects of aerosol emissions during the 21st century. However, we also conducted a historical simulation from year 1850–2000. First, we simulated the aerosol forcing with ECHAM-HAMMOZ at 10 year intervals. For each of these time points, we ran ECHAM-HAMMOZ continuously for three years with an additional three-month spin-up period. This approach is computationally cheaper than running the model for each year with different emissions and also reduces noise from inter-annual variability in the results. Then, we calculated the mean effective radiative forcing (ERF) including both shortwave and longwave contributions for each calendar month from the three-year runs as the mean difference in net radiation at the top-of-the-atmosphere between a given emission scenario and a reference run with year-1850 anthropogenic and non-agriculture biomass burning aerosol emissions (ERF takes into account fast feedbacks such as changes in cloud cover and is thus a better predictor of temperature change (Haywood et al 2009, Lohmann et al 2010, Myhre et al 2013b)). In this way, ERF for each month and grid-cell for a given year was calculated as an average of the same month over three consecutive years. Due to the interpolation, the temporal patterns showed less year-to-year variability than it would have been with every year simulated individually. Mineral dust and sea salt emissions were calculated online identically in all simulations. We used linear interpolation to calculate the aerosol forcing for the years between the 10 year intervals, and applied that forcing into the one-layer atmosphere of the UVic ESCM from 1850 until 2000; the same 10 year linear interpolation approach was also used for simulations covering the 21st century.

For simulations over the years 2001–2100 with the UVic ESCM, we chose RCP4.5 (Thomson et al 2011) as the baseline scenario. Thus, emissions of CO₂ and forcings from non-CO₂ greenhouse gases came from the RCP4.5 scenario in all of our UVic ESCM simulations and the only difference between the runs was aerosol forcing. To explore the range of different aerosol emissions across RCPs, we did simulations using aerosol forcing (calculated with ECHAM-HAMMOZ) from RCP2.6 (van Vuuren et al 2011b), RCP4.5, or RCP8.5 (Riahi et al 2011). We designated these scenarios as 2.6AER, 4.5AER, and 8.5AER, respectively. We excluded RCP6.0 as its emissions fall between those of RCP4.5 and RCP8.5, and we only wanted to cover the range of possible aerosol emissions in the RCPs and not necessarily study each individual RCP.

Since the RCPs do not span the full range of conceivable aerosol emission scenarios (Rogelj et al 2014a, Fiore et al 2015), we built two idealized sensitivity scenarios to represent a very aggressive aerosol mitigation scenario and a scenario where aerosol emissions are allowed to grow with increasing CO₂ emissions. These new sensitivity scenarios allow us to assess climate and health effects from a wider range of possible outcomes. The sensitivity scenarios are based on RCP4.5, modifying only the anthropogenic aerosol emissions (i.e. keeping the same emissions from biomass burning, sea salt, and mineral dust as in 4.5AER). For each pair of sector and species of anthropogenic aerosol emissions (e.g. BC from agricultural waste burning or SO₂ from land transport), we calculated the ratio of the aerosol emissions and global CO₂ emissions at each year. For the high-end sensitivity scenario (designated as HIGH), we assumed that these ratios stayed at the level of the year 2005 throughout the 21st century. For the low-end sensitivity scenario (designated as LOW), we assumed that these ratios would reach a minimum level in 2030 (i.e. minimum ratio of aerosols per global CO₂ emissions found across RCPs 4.5, 6.0, and 8.5) and stay at that level afterwards. We did not consider the ratios in RCP2.6 as the global net CO₂ emissions become negative in this scenario. The ratios for the years 2010 and 2020 were linearly interpolated using values of 2005 and 2030. To calculate the actual aerosol emissions for HIGH and LOW in a given year, we multiplied the previously calculated ratios with global CO₂ emissions for that year and then distributed these emissions spatially based on RCP4.5 for that same year. Figure 1 shows the global total anthropogenic emissions of SO₂, BC, and OC for all aerosol emission scenarios used in this study. Scaling aerosol emissions with CO₂ emissions is based on assumptions that aerosols are co-emitted with CO₂ and that sectoral CO₂ emissions follow the global trend. Even though these assumptions are less valid for some sectors, LOW and HIGH scenarios appear appropriate in the context of the current global-level sensitivity study of different aerosol emission pathways.

Zoom In Zoom Out Reset image size

We used the surface aerosol concentrations from ECHAM-HAMMOZ simulations to calculate premature mortality caused by chronic exposure to ambient PM_2.5 resulting from anthropogenic aerosol emissions. We mainly followed the methodology by Ostro (2004) with cause-specific relative risk calculated with the Integrated Exposure-Response (IER) model by Burnett at al (2014). For the reference (counterfactual) scenario, we made a simulation without even the minimal year-1850 anthropogenic aerosol emissions, but with year-1850 biomass burning, sea salt, and mineral dust aerosol emissions (this differs from the reference scenario for ERF calculations, which had year-1850 anthropogenic and non-agriculture biomass burning aerosol emissions as explained in section 2.2). For each simulation, we diagnosed the three-year mean surface PM_2.5.

We estimated the PM_2.5-caused mortality for ischemic heart disease (IHD), cerebrovascular disease (stroke), chronic obstructive pulmonary disease (COPD), and lung cancer (LC). To get the anthropogenic contribution, we subtracted the mortality calculated from the counterfactual scenario from the results of other scenarios. We calculated the total annual premature mortality rate E (deaths per year per unit area) using the following formula:

$$\begin{equation} {\rm{E}} = \mathop \sum \nolimits_{{{i}},{{j}}} {\rm{A}}{{\rm{F}}_{{{i}},{{j}}}} \times {{\rm{B}}_{{{i}},{{j}}}} \times {{\rm{P}}_{{i}}} \end{equation} \tag{ 1 }$$

where AF_i,j is the attributable fraction in age group i for cause j, B_i,j is the baseline mortality rate (deaths per year per capita), and P_i is the exposed population in age group i. AF_i,j is 1−1/RR, where relative risk is based on the IER model by Burnett et al (2014). Burnett et al (2014) base their model on the following equation:

$$\begin{equation} {\rm RR}_{i,j} (z){\kern-2pt} ={\kern-2pt} 1+\alpha_{i,j} \left\{{\kern-2pt} 1-{\rm exp} \left[-\gamma_{i,j} (z-z_{{\rm cf},_{i,j}})^{\delta_{i,j}}\right]\right\} \end{equation} \tag{ 2 }$$

where z is the PM_2.5 concentration (in μg m⁻³) in a given year in a given aerosol emission scenario and z_cf is the counterfactual concentration that is the lower limit for any adverse health effects. Burnett et al (2014) provide an ensemble of 1000 sets of parameters α, γ, δ, and z_cf for each affected age group i and cause j. Following the example by Apte et al (2015), we evaluated equation (2) with all different ensemble values at PM_2.5 concentrations up to 409.9 μg m⁻³ with an interval of 0.1 μg m⁻³, and created a lookup table for the ensemble mean and upper and lower bounds of the 95% confidence interval (defined by 2.5% and 97.5% quantiles) of RR. The use of ensemble mean and lower and upper bounds relies on the assumption that the 1000 parameter sets provided by Burnett et al (2014) are equally probable.

Projected changes in demographics will have a major impact on PM_2.5-caused mortality, even if the air quality itself does not change (Silva et al 2016). Since we wanted to focus on the health implications of changing aerosol concentrations only, we used constant baseline mortality rates in six different world regions (WHO 2008) (see figure 2 of Partanen et al (2013) and year-2010 population density data (SEDAC 2005) for all time intervals. Using the same baseline mortality rate for all countries within a single WHO region hides some spatial variation, but on a global level the differences partly cancel each other and we consider these limitations to be acceptable in the context of this global-scale sensitivity study. WHO (2008) provides baseline mortality rates and fractions of age groups of 30–44, 45–59, 60–69, 70–79, and 80+. We calculated the mortality caused by IHD and stroke separately for these age groups and by COPD and LC for population over 30 years of age. Since Burnett et al (2014) provide relative risks for more fine-grained age groups, we used the average of their relative risks to calculate the relative risk for each of our age group (e.g. relative risk for 30–44 was calculated as the average of age groups of 30–34, 35–39, and 40–44).

Zoom In Zoom Out Reset image size

Figure 2 shows the global mean ERF from changes in aerosol emissions through time (although the results include also biomass burning emissions, most of the effect comes from changes in anthropogenic aerosols). The maximum negative ERF of −1.41 W m⁻² in the historic run occurred in 1990. In 2011, we obtained values between −1.19 and −1.30 W m⁻² for the three AER simulations, which compares to the best estimate of −0.9 (−1.9 to −0.1) W m⁻² from the latest IPCC report for the aerosol ERF in 2011 relative to 1750 (Myhre et al 2013b). The year-2000 global mean ERF (−1.31 W m⁻²) was also somewhat stronger than the multi-model mean of −1.17 W m⁻² obtained by Shindell et al (2013). To assess the effect of interannual variability on ECHAM-HAMMOZ results, we calculated the annual ERF for each of the three individual years used to compute the mean year-2000 ERF (as explained in section 2.2, results for each 10 year interval were based on running ECHAM-HAMMOZ continuously for three years). The annual ERF ranged between −1.28 and −1.36 W m⁻². Thus, the effect of ECHAM-HAMMOZ interannual variability appears substantially smaller than inter-model uncertainty in aerosol ERF.

The transient global mean ERF was fairly similar for the three AER simulations. The magnitude of aerosol ERF decreased steadily to between −0.15 W m⁻² and −0.39 W m⁻² in 2100. The largest difference across the AER simulations was 0.34 W m⁻² in 2040. These results are similar to results by Westervelt et al (2015). The differences were much larger for our low- and high-end estimates. In LOW, the magnitude of aerosol ERF decreased quickly and was only −0.14 W m⁻² in 2030. The subsequent decline (driven by decreasing global CO₂ emissions) was slower, reaching −0.02 W m⁻² by 2100. In HIGH, on the other hand, the magnitude of aerosol ERF continued to increase until 2040 due to the increase of CO₂ emissions, reaching a peak value of −1.63 W m⁻². After 2040, the magnitude of aerosol ERF started to decline but it was still −0.82 W m⁻² in 2100.

Aerosol forcing was strongest in the northern hemisphere and especially over Eastern Asia and over oceans nearby the shipping routes. In 2005 for 4.5AER, the aerosol forcing was strongest in the northern hemisphere and especially over Eastern Asia and over oceans, reaching values of about −12 W m⁻² in these regions (figure S1(a)). In 2100 for 4.5AER, northern hemisphere oceans, East Asia, central Africa and the South Atlantic Ocean near the African coast had noticeable negative aerosol ERF (figure S1(b)). In the same year for LOW, most areas had negligible ERF as negative and positive values mostly took place next to each other (figure S1(c)). However, the east coast of North America experienced positive aerosol ERF in 2100 for both 4.5AER and LOW. This positive forcing was most likely caused by the substantial reduction in biomass burning emissions over the region. The forcing pattern in 2100 for HIGH (figure S1(d)) was generally similar to, although larger in magnitude than, 4.5AER. The spatial patterns of forcing in 2.6AER and 8.5AER were also similar to 4.5AER.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As with aerosol forcing, the effects on global mean surface temperatures across AER scenarios were similar, but results differed considerably for LOW and HIGH (figure 3(a)). The global mean surface temperature difference between the AER simulations remained less than 0.18 °C, year-2100 difference being 0.17 °C. The difference in end-of-the-century (2096–2100) global mean surface temperature between 2.6AER and 8.5AER was 0.17 °C, comparable to the corresponding aerosol-caused temperature differences between RCP2.6, RCP4.5, and RCP8.5 scenarios by Westervelt et al (2015) (0.07 °C–0.23 °C). On the other hand, the difference between LOW and HIGH was 0.86 °C at its maximum (in 2061) and 0.61 °C in 2100. It is also notable that simulated warming diverges quite rapidly over the years 2015–2040 in the LOW and HIGH scenarios, which is not the case among the AER scenarios. Previous analyses of contributions to uncertainty in future warming projections have suggested the scenario uncertainty remains quite low over the coming decades (e.g. Booth et al 2013); by contrast, our results suggest that considering a wider range of plausible aerosol scenarios would considerably increase the contribution of scenario uncertainty to near-term temperature projections.

We also calculated global mean warming rate as a 10-year running mean of global mean surface temperature change over time (figure 3(b)). In all simulations, there was an initial peak in the warming rate (about 0.4 °C per decade) in 2002 due to rapid warming in the last years of the historical simulation (not visible in figure 3(a)). Among the AER simulations, 2.6AER had the highest warming rate with a maximum value of 0.34 °C per decade in 2037. For most of the second half of the century, 4.5AER had the highest warming rate among AER scenarios. In LOW, the warming rate increased quickly during the first decades with rapidly decreasing aerosol emissions, reaching a maximum value of 0.48 °C per decade in 2032. This was 81% higher than the warming rate with 4.5AER in the same year. On the other hand, in HIGH, the increase of aerosol negative ERF until 2040 and the gradual decrease afterwards (figure 2) resulted in the most stable warming rate of all the simulations, remaining below 0.2 °C per decade between 2008 and 2058.

Figures 3(c) and (d) show the spatial differences in mean surface temperature for 2090–2100 between LOW and 4.5AER, and between HIGH and 4.5AER, respectively. These differences were very similar with respect to magnitude and spatial patterns, but had opposite signs. The largest differences were in high-latitude regions, particularly in the northern hemisphere, with surface temperature difference between LOW and HIGH reaching 1.3 °C at the North Pole. Although ERF was highly localized in specific regions, the surface temperature differences created by different aerosol emissions were fairly smooth, particularly longitudinally.

We calculated the anthropogenic surface PM_2.5 contribution as the difference between a given scenario and the preindustrial year-1850 simulation without any anthropogenic emissions. The largest anthropogenic PM_2.5 surface concentrations in year 2000 were in China, where the concentrations were between 25 and 50 μg m⁻³ (figure S2(a)). The anthropogenic PM_2.5 concentrations in year 2000 were between 5 and 25 μg m⁻³ over large parts of North America, Europe, Asia, and Northern Africa. To compare our results with previous studies, we calculated also the difference between years 2000 and 1850 including the anthropogenic aerosol emissions in both years (figure S2(b)), and the year-2005 surface PM_2.5 concentrations for the 4.5AER simulation (figure S2(c)). When compared to the model results of Silva et al (2013) (their figure S3), our simulated increase in surface PM_2.5 concentration between 1850 and 2000 (figure S2(b)) was similar or slightly larger in North America and Europe, comparable in South America and the polluted regions of Asia, slightly higher in Northern Africa, and slightly lower in parts of Central Africa. PM_2.5 concentrations in 2005 for 4.5AER (figure S2(c)) were largely similar to those of Brauer et al (2012) (their figure 2), who combined modeling to ground and satellite measurements. The largest differences were in Northern Africa and Australia, where our model predicted higher PM_2.5 concentrations, most probably due to differences in dust emissions.

Figure 4(a) shows global total premature mortality due to PM_2.5 in the different aerosol emissions scenarios for the central estimates of the RR values (see below for an uncertainty analysis). As discussed in section 2.3, these numbers represent what the mortality would be with present-day population density and baseline mortality rates, and cannot therefore be directly interpreted as future predictions. Taking into account the combined effect of population growth and changing baseline mortality rates would give higher global total future premature mortality (Silva et al 2016). Larger mortality in all simulations would also mean that the differences between our scenarios would be larger. For example, LOW would lead to larger avoided mortality when compared to the other scenarios. Geographically, a larger fraction of global total mortality would take place in Africa due to the fastest projected population growth.

The PM_2.5-caused mortality for 4.5AER in 2010 was about 2 371 800 (95% confidence interval of 1 329 300–2 925 200) deaths per year. As aerosol emissions decreased in all AER scenarios, premature mortality due to PM_2.5 also decreased steadily and similarly across scenarios, and was only about 525 700 (258 400–757 300) deaths per year for 4.5AER in 2100. 2.6AER had the lowest mortality across all AER scenarios. In HIGH, mortality increased until 2030, when it reached 2 780 800 (1 607 000–3 366 300) deaths per year. HIGH reached its minimum mortality in 2080 (1 150 900 (608 800–1 530 800) deaths per year), followed by a small increase as the emissions increased slightly until the end of the century. PM_2.5-caused premature mortality in LOW decreased quickly along with aerosol emissions. It was below 300 000 deaths per year from 2030 onwards and went very close to zero (below 30 000 deaths per year) from 2080 onwards. Throughout the simulations, the uncertainty bounds of the AER simulations were overlapping with each other (table S1). On the other hand, the higher bound of LOW was below the lower bound of HIGH from year 2030 onwards. This shows that the AER scenarios are fairly similar to each other, whereas LOW and HIGH differ from each other significantly.

Even though the global mortality in all simulations was positive throughout the 21st century, the mortality was below zero in some regions near the end of the century. For example, in LOW, the mortality in 2100 was negative in India, eastern parts the United States, and in parts of South America (figure S3(a)). The negative (i.e. avoided) mortality means that in LOW, PM_2.5 concentrations in these regions were below the reference level from year-1850 without anthropogenic aerosol emissions. This happened because biomass burning aerosol emissions in LOW, which were prescribed according to RCP4.5, were in some regions lower than in the year-1850 reference simulation. The largest difference in PM_2.5 was due to differences in organic carbon emissions (figure S3(b)).

Looking at the spatial pattern of mortality, the differences between LOW (figure 4(b)) or HIGH (figure 4(c)) vs. 4.5AER were mostly mirror images of each other; regions with a negative difference in mortality for LOW had a positive difference of similar magnitude for HIGH. The differences were largest in heavily populated areas of Asia and Europe.

Our estimate of the present-day (2000) mortality (2 263 500 (1 267 500–2 815 300) deaths per year) was about 8% higher than the estimate of 2.1 million deaths per year by Silva et al (2013) (years 2000 vs. 1850 in their study) and 33% higher than the estimate of 1.7 million deaths per year by Silva et al (2016). The main reason for the higher mortality in our study was most likely that we defined the reference scenario completely without anthropogenic aerosol emissions, whereas Silva et al (2013, 2016) used year-1850 emissions that include a small anthropogenic contribution. Therefore, our estimate for anthropogenic PM_2.5 concentration was higher over North America and Europe (figure S2(a)) than in Silva et al (2013). To better compare with the results by Silva et al (2013, 2016), we recalculated the mortality by including anthropogenic emissions in the year-1850 reference simulation. This resulted in mortality of 1 831 700 (1 067 300–2 178 900) deaths per year for the year 2000, which is very close to the estimate by Silva et al (2016), who used the same concentration-response function as we did. On the other hand, our present-day PM_2.5-induced mortality estimate is smaller than the 3.15 (1.52–4.6) million deaths per year estimated by Lelieveld et al (2015). Our results can also be compared to those of Likhvar et al (2015) who estimated that maximum feasible reductions of aerosol emissions would lead to 1.5 (0.4–2.4) million avoided deaths per year in 2030 compared to 2010. A comparable figure in our study is the difference in mortality between 4.5AER in 2010 and LOW in 2030 (2.1 (1.2–2.5) millions of deaths per year).

Johnston et al (2012) found that results varied by a factor of 2.3 when using two different functions to quantify PM_2.5-induced mortality resulting from landscape fires. To test the sensitivity of our results to the choice of the concentration-response function, we calculated the relative risk also using the functions given by Ostro (2004) (see supplemental material available at stacks.iop.org/ERL/13/024028/mmedia). These results did not differ considerably from the results using the IER model by Burnett et al (2014) (compare tables S1 and S2). When using Ostro (2004) to calculate relative risk, the annual global total mortality was between 12% and 27% higher than with the IER model; also, 95% confidence intervals were wider with Ostro (2004) than with the IER model.

Our premature mortality results are likely to be biased low due to the coarse model resolution. Kodros et al (2016) showed that decreasing model resolution from 2° × 2.5° to 4° × 5° decreased by 16% global total mortality caused by aerosols from domestic combustion. Lower mortality with coarser resolution is caused by potential dilution of PM_2.5 in larger grid cells as well as reduced overlap of densely populated areas and high PM_2.5 concentrations (Kodros et al 2016).

As discussed above, ambitious reductions in aerosol emissions resulted in both increased warming and greatly decreased PM_2.5-induced premature mortality. This potential trade-off is highlighted in figure 5, where we plot the trajectories of surface temperature change and premature mortality. Although the surface temperature increase and premature mortality decrease were faster in LOW, the situation in year 2100 was not very different from that in the AER simulations. On the other hand, the trajectory in HIGH was considerably different, especially with increasing premature mortality during the first decades of the 21st century.

Zoom In Zoom Out Reset image size

We assessed the premature mortality from changes in aerosol emissions only, omitting possible feedbacks that climate change could have on air quality (Fiore et al 2015, Silva et al 2013, von Schneidemesser et al 2015, Allen et al 2016). Greenhouse gas emission reductions can also lead to decreases in aerosol emissions (West et al 2013), which mitigates the trade-offs presented here. Aerosol emissions could also affect human mortality indirectly; in particular, heat waves are known to increase premature mortality (Mitchell et al 2016) and aerosol-induced cooling could reduce heat wave intensities (Péré et al 2011). Comparing premature mortality due to heat waves to premature mortality due to PM_2.5 exposure would prove particularly interesting and should be the topic of future studies. Unfortunately, existing exposure-response functions for mortality caused by heat waves require city-specific coefficients and daily mortality data that are not available for a global study. For example, there is evidence that the temperature threshold above which heat waves cause mortality can vary by more than 10 °C across specific cities (Baccini et al 2008), so extrapolating the relationships available for a handful of cities to the global scale appears risky at best. Moreover, the coarse spatial and minimal internal variability of the UVic ESCM would make it difficult to diagnose heat waves from the model results. Simplistic back-of-the-envelope calculations could be made at a global level by assuming linear relationship between global mean surface temperature change and climate-change caused mortality (Löndahl et al 2010), but we deemed this approach too unreliable to calculate meaningful numbers.