A number-based inventory of size-resolved black carbon particle emissions by global civil aviation

With the rapidly growing global air traffic, the impacts of the black carbon (BC) in the aviation exhaust on climate, environment and public health are likely rising. The particle number and size distribution are crucial metrics for toxicological analysis and aerosol-cloud interactions. Here, a size-resolved BC particle number emission inventory was developed for the global civil aviation. The BC particle number emission is approximately (10.9 ± 2.1) × 1025 per year with an average emission index of (6.06 ± 1.18) × 1014 per kg of burned fuel, which is about 1.3% of the total ground anthropogenic emissions, and 3.6% of the road transport emission. The global aviation emitted BC particles follow a lognormal distribution with a geometric mean diameter (GMD) of 31.99 ± 0.8 nm and a geometric standard deviation (GSD) of 1.85 ± 0.016. The variabilities of GMDs and GSDs for all flights are about 4.8 and 0.08 nm, respectively. The inventory provides new data for assessing the aviation impacts.


Supplementary Note 1. Data preparation
A database has been prepared as shown in Supplementary Table 1 to investigate the correlation between the GMD/GSD of BC particles and the engine parameters. Most of the data are measured by utilizing the standardized system, e.g. the Missouri University of Science and Technology (Missouri

Supplementary Note 3. GMD-T 3 correlation
The abundant ground measurements can be utilized to estimate the particle size at cruise, if the In the test experiments, the BC particle size distributions were measured for a typical turbofan engine operated at different simulated altitudes (from sea-level-static to 15.2 km, about 0.8 Mach) in a test facility at Arnold Engineering Development Center (AEDC) 63 . The dataset showed that the size distribution of BC particles was strongly dependent on the engine inlet temperature T 3 . The fitted relation between the particle size and T 3 is applicable for different simulated altitudes inside the test facility. There is no significant correlation between the residuals and the altitudes.
The analysis of variance (ANOVA) test was conducted to investigate the differences among the residuals at different altitudes, which indicates that there is no statistically significant difference among them . The ANOVA table is shown in Supplementary Table 10, with a p value of 0.73, which means there is a high risk (73%) to claim that the residuals are different. As a result, the experimental results indicate that the dependence of BC particle size on T 3 is applicable to different simulated altitudes.
In the present study, the T 3 temperature is estimated following the method in FOX method 26 , which is briefly described in Supplementary Note 18. Supplementary Figure 1(b) shows the relation between GMDs and estimated T 3 temperatures. There is no distinct difference between the cruise and ground data, therefore a quadratic equation (S.2)   41   2  3  3  1  2  3 1000 1000 is adopted to fit the data. All the data are within the prediction interval of 95% confidence. The least square linear regression was utilized to estimate the parameters.
In Supplementary Figure 16, the data at different simulated altitudes were compared with the data utilized to develop the GMD-T 3 correlation. The results indicate that new dataset (colorful dots) agrees well with the utilized dataset (grey dots). All of the new data are within 95% prediction interval of the fitted correlation.

Supplementary Note 4. GMD-EI m correlation
The mass emission index EI m is also a potential variable of a universal fitting relation for the GMDs The least square linear regression was utilized to estimate the parameters. Two statistical metrics, adjusted R 2 , and root mean square deviation (RMSE), are utilized to evaluate goodness of the fitting results. The adjusted R 2 was utilized the influence of the number of parameters in the fitted relation.

Supplementary Note 5. GMD-T 3 &EI m correlation
The performances of the above three different fitted relations are compared in the bottom row of Supplementary Figure 1. The often adopted GMD-fuel flow rate (thrust) relation has the lowest quality, with a high RMSE (4.03) and a relatively low adjusted R 2 (0.63). The RMSEs respectively decreased by 8.8% and 30.0% when using the GMD-T 3 and GMD-EI m relations. The adjusted R 2 are also improved by 7.9% and 28.6%. 42 Although the qualities of the fitted correlations have been improved based on T 3 and EI m , the data scattering is still significant, which indicates that single type input, such as the fuel flow rate, thrust, or T 3 is insufficient to adequately predict the GMD. Therefore, we try to develop the correlation for The improvement of the new correlation is discussed in the main article.
The coefficients of all the fitted relations are shown in Supplementary Table 3. As shown in the main body, the quality of the GMD-T 3 &EI m correlation outperforms the others, as evidenced by the higher 43 adjusted R 2 and the lower RMSE. The bias-variance trade-off tests were also conducted to investigate whether the GMD-T 3 &EI m correlation was over-fitted in Supplementary Note 6.

Supplementary Note 6. Bias-variance trade-off tests
In order to investigate the bias-variance trade-off of our models, all the data of measured particle size distributions were randomly divided into the training (33 data points about 60%) and validation (20 data points about 40%) sets. One thousand stochastic groupings were generated and calculated for the bias-variance trade-off tests to avoid the influence of the grouping method. Both the training and validation datasets contained the data from all the engines.
Eight different models were utilized as shown in Supplementary Table 4 Through the bias-variance trade-off tests and comparison with the independent APEX dataset, we demonstrate that our model does not over-fit the data.

Supplementary Note 8. GSD estimation
The Geometric Standard Deviation (GSD) of the BC particles is also an important parameter for the particle number emission estimation. The relations between GSD and EI m , T 3 are fitted. Different from the situation for GMD, the dependence of GSD on T 3 is not strong enough, as shown in Supplementary Figure 4 (a). The data point has large scattering around the fitted line. The correlation between GSD and EI m is much stronger. A simple linear relation between the GSD and the EI m is adopted to predict the GSD in this study, as shown in Supplementary Figure 4 (b). The variation of GSD is not as large as that of GMD. As a result, the fitted GSD-EI m correlation is adopted in this study.

Supplementary Note 9. The uncertainties of particle mass parameters
The uncertainties of ε shown in Abegglen et al. 39 were the 95% confidence interval, not the standard deviation. Besides, the uncertainties were significantly overestimated. In Abegglen et al. 39 , ε was estimated by fitting the power-law relationship between particle mass and size pm m C d   , (S.5) The particle mass was directly measured by Centrifugal Particle Mass Analyzer (CPMA; 46 Cambustion Ltd., Cambridge, UK) at fixed particle sizes and engine thrusts. Based on the measurements, the standard deviations for the mass of the size-and thrust-selected particles were maximum 16.6% for the particles smaller than 40 nm and 11.0% for those larger than 40 nm. We re-analyzed the data in the reference to determine the uncertainties of the mass data more accurately.
Supplementary Figure 5 shows the measured particle mass and uncertainties (error bars) for 33%, 67% and 105% thrusts from Abegglen et al. 39  The shaded area is the results estimated using an uncertainty of 0.025 for ε. The estimated uncertainties are close to those of the measurements.
According to Abegglen et al. 39 , the variation of prefactor C is very small, within about 5%. As a result, for a given particle size (d m ), the change of mass caused by the parameter uncertainties is where E(m p ) and E(ε) are the means of particle mass and ε. For the particles with an electrical mobility diameter (d m ) between 20 to 150 nm, the changes of particle mass are about ±40% with the ε uncertainty utilized in this study (∆ε=0.025). These estimated particle mass uncertainties are already larger than the measurements (11.0% to 16.6%). We assume that it is reasonable to use larger uncertainties for particle mass considering the measurements are only for one type of engine.
However, the changes of particle mass can be as large as -90% to 2000% if the uncertainties of ε

Supplementary Note 10. Residuals of the EI n estimation
Supplementary Figure 6 shows the residuals of all the Monte Carlo runs separately for the entire data (cruise + ground), cruise and ground measurement data. The residuals are close to normal distributions. The mean of the residual for all the data is 0.84×10 13 /kg-fuel with a 95% confidence interval between (−0.16 1.84)×10 13 . The mean residual is small compared to the BC particle number emission index which is normally at the magnitude of 10 14 /kg-fuel. The standard deviation is 1.43×10 14 /kg-fuel, the same magnitude with the emission index.
The EI n data for the cruise and ground conditions are respectively over-and underestimated, evidenced by the mean residuals of 7.82×10 13 and −1.39×10 13 /kg-fuel. The bias may be caused by both the experimental and modeling uncertainties, e.g. the flight parameters, atmospheric conditions and measurement technologies. The limited number of data points may also influence the results.
The biases are within the uncertainties of the global mean EI n (BC) estimations, which have the standard deviations of 1.00×10 14 /kg-fuel for LTO and 1.22×10 14 /kg-fuel for CCD as shown in Supplementary Table 6.

Supplementary Note 11. Sensitivity tests for the influences of BC mass emissions
The BC mass emission directly influences the estimation of particle number emission. Sensitivity tests were conducted to investigate the influence of BC mass emissions on the particle number emission and the size distribution. The estimations based on the EEA dataset were utilized as the baseline in the sensitivity tests. Besides the baseline, four different cases with substantial changes in BC mass emissions were included. It was assumed that the BC mass emissions were respectively 33% (1/3, 3.2 Gg), 50% (1/2, 4.8 Gg), 150% (14.3 Gg) and 200% (19.1 Gg) of the baseline emission (9.5 Gg). The range of these values is able to cover most of the BC mass emissions in the literature, as shown in Table 2. 48 Supplementary Figure 9 shows the results of the sensitivity tests. The changes of particle number emissions ( Supplementary Figure 9 a) are much smaller than the corresponding changes of mass emissions, because higher EI m (BC) induces larger and heavier particles, which alleviates the increase of BC particle number and makes the particle number emission relatively stable. It is similar when EI m (BC) is lower with smaller and lighter particles. Supplementary Figure 9 a) shows that, for the cases with 50% and 150% mass emissions, the resulted particle number emissions are respectively In summary, the sensitivity tests indicate that our estimations of the BC number emission and size distribution were robust, even when we applied substantial changes of the BC mass emission from 33% to 200% of the baseline. The sensitivity tests demonstrate that the BC particle number emission and size distribution are rather robust, considering the substantial changes of the BC mass emissions used in the tests. Therefore we consider the uncertainties are at reasonable based on the results.

Supplementary Note 12. Monte Carlo simulations and uncertainty estimations
By The mean and standard deviation from the MC runs are shown in Figure 4. The GMD and GSD of the global BC size distribution were obtained by fitting the particle numbers in the size bins into the lognormal distribution for each MC run, and the uncertainties of GMD and GSD were estimated from these MC runs.
In Supplementary Table 8, the variabilities were calculated as the standard deviations of the emissions from all the 27 million flights at the seven sub-phases (taxi, take-off, climb-out, approach, climb, cruise and descent). The uncertainties are about 17% to 21% of the variabilities.

Supplementary Note 13. Estimation of BC particle number emission from open burning
The The mass to number conversion follows the similar methods for the aviation emissions. Assuming that the particle is spherical, we can express the particle mass m p as where ρ is the effective density of the BC particle, which is assumed to be 10 3 kg m -3 ; d is the volume equivalent diameter, which follows the lognormal distribution as shown in Supplementary   Table 9), which was mainly due to the inclusion of emissions from kerosene wick lamps, gas flaring, and re-estimation of the emissions from China using the regional coal statistics 56 .
. (S.14) The prefactor contains the density of the primary particle. However, Abegglen et al. 39 did not assume any density. They directly measured the particle mass with desired sizes selected by a Differential Mobility Analyzer (DMA) under specific thrusts. The ε and prefactor C were estimated by fitting the power-law relationship between particle mass and size. It is difficult to estimate ρ 0 from the prefactor C, because the parameters k and f were not specifically estimated.

Supplementary Note 15. Particle mass-number conversion
The particle mass emission index can be expressed as the integration of all the mass of the single  The new data (red dots) agree well with the utilized dataset (grey dots). All of the new data are within 95% prediction interval of the fitted correlation. At 31% thrust, the main annulus started, leading to lower BC mass emission. The measured GSD at 31% is smaller than that at 17%, which indicates a strong relation between GSD and EI m (BC).
In general, the proposed method successfully captures the unique changes of GMD of the BC particles generated in this double annular engine, even though the method is developed based on the engines with a single annular combustor. The results demonstrate the necessities to include both the T 3 and EI m (BC) as the predictors for GMDs. They also show the potential abilities of the proposed method to estimate the BC particle size distribution and the number emission for significantly different engine designs and operating conditions.

Supplementary Note 17. R 2 for the nonlinear regression of GMD-T 3 &EI m correlation
Normally, the adjusted R 2 is applicable to the linear fittings. Here we show that it is also valid for the current GMD-T 3 &EI m correlation. In order to calculate an effective  For the linear regressions with a constant term, the two assumptions can be satisfied. Equation (6) meets the first assumption, because the Levenberg-Marquardt nonlinear least squares algorithm was utilized. For the second assumption, scaling (multiplied by a) of the model f cannot improve the fitting, and shifting (added by b) of the model f also barely changes the results, because there is a constant term A 4 · A 6 . As a result, R 2 is still applicable for the nonlinear regression of GMD-T 3 &EI m correlation.

Supplementary Note 18. T 3 estimation
The T 3 temperature is estimated following FOX method 26 . The temperature at the ground condition can be estimated as: where P 3 is the combustor inlet pressure with the unit of atm, 00  is the overall pressure ratio, 2 T and 2 P are respectively the engine inlet temperature and pressure, namely the ambient temperature and pressure during the ground experiments,  is the heat capacity ratio (  = 1.4) and P  is the polytropic efficiency which is assumed as 0.9 26 .
For the cruise condition, the engine inlet temperature and pressure have to be revised based on the speed of aircraft (Mach number, Ma) and the ambient pressure and temperature at cruise (P s , T s ): The EEA data provides the total emissions respectively for the whole LTO and CCD phases. The flight conditions (e.g. fuel flow rate, speed and altitude) are not constant during LTO or CCD. The different flight conditions influence the BC emission and mass to number conversion. In order to better estimate the number emissions, the EEA mass emission data for LTO and CCD phases are further divided into seven sub-phases, taxi/take-off/climb-out/approach for the LTO cycle and climb/cruise/descent for the CCD phase.
In order to separate the data, we first estimate the durations of the sub-phases at various flight