Four-dimensional aircraft emission inventory dataset of the landing-and-takeoff cycle in China (2019–2023)

The rapid growth of the aviation industry has resulted in aircraft emissions during landing and takeoff (LTO), which have direct and increasingly adverse impacts on air quality and human health. An accurate and high-resolution LTO emission inventory is crucial for investigating these adverse effects, with the LTO emission having unique three-dimensional (3D) spatial characteristics and typical hourly temporal variations. This study integrated the emission calculation and flight trajectory recognition methods to establish a four-dimensional (4D) aircraft emission inventory dataset of China's LTO cycle (4D-LTO emission inventory dataset) from 2019 to 2023. The dataset has a high spatial–temporal resolution (hourly, 0.03° × 0.03° × 34 height layers) and incorporates calculation emissions accurately. Moreover, the actual taxi-out/taxi-in time for each flight was determined by a statistical model of taxi time and some aircraft in schedule based on 38 million flights. Each flight's climb/approach time was also obtained based on mixing layer height (MLH) and the height–time nonlinear relationship. Additionally, we calculated the LTO emission for China's flight, establishing the hourly emission inventory based on each mode's running time, emission index, and fuel flow. We obtained the flight trajectory core of each airport based on measured flight trajectories and the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) to depict the spatial distribution. Then, each flight's takeoff/landing direction and trajectory were identified from the wind direction and relative departure/arrival airport position. The findings indicate that during COVID-19, the LTO number in 2020–2022 was reduced to 73.1 %, 77.6 %, and 48.7 % of 2019 levels, respectively. However, by 2023, the LTO number has rapidly bounced back to 95.3 % of 2019 levels. The recovery rate during daytime (06:00–23:00 UTC+8) was 41.6 % higher than during nighttime (00:00–05:00 UTC+8). The emissions of various pollutants were measured as follows: hydrocarbon (HC), carbon monoxide (CO), nitrogen oxides (NO_x), particular matter (PM), and sulfur dioxide (SO₂) are 3.2, 46.1, 62.3, 1.1, and 18.4 Gg. LTO emissions' horizontal characteristic is the distance along the runway and spread. This elongated distribution will be hidden if a rough grid is used (e.g., 0.36° × 0.36°) and the emissions are evenly distributed. Moreover, LTO emission height characteristic decreases with height, and the maximum height varies with MLH. Emissions above the standard height set by the International Civil Aviation Organization standard height (∼ 915 m) are not estimated. For example, NO_x emissions above 915 m during various months make up an average of 24.6 % (9.9 %–37.5 %) in the LTO cycle, indicating the emissions are significantly underestimated when using the ICAO method. Compared with conventional spatial allocation methods, our dataset provides a more accurate representation of the actual LTO situation in both the horizontal direction and height at different times. Our 4D-LTO emission inventory dataset and its adaptable methodology are valuable resources for researching temporal and spatial variations, air quality, and health impacts of aircraft emissions in the LTO cycle. The dataset can be accessed via https://doi.org/10.5281/zenodo.13908440 (Lang et al., 2024).

1 Introduction

The aviation industry has experienced rapid growth in recent years. However, aircraft emit pollutants such as NO_x, CO, SO₂, HC, and PM during operation; affect air quality; and have adverse effects on human health and human life (Wang et al., 2022; Dissanayaka et al., 2023; Pandey et al., 2024). It has been estimated that 8000–58 000 premature mortalities each year are attributable to aviation emissions (Barrett et al., 2010; Eastham and Barrett, 2016; Quadros et al., 2020; Eastham et al., 2024). Establishing an accurate aircraft pollutant emission inventory is crucial to investigating the impact of aircraft emissions on the environment and health.

According to the standard height (∼ 915 m) of the mixed layer height (MLH), the International Civil Aviation Organization (ICAO) divides the flight process of the aircraft into the landing-and-takeoff (LTO) cycle phase and climb–cruise–descend phase (Kurniawan and Khardi, 2011; Bao et al., 2023). The LTO cycle occurs near the ground and affects the air quality near the airport and the health of the surrounding residents (Christodoulakis et al., 2022). Therefore, many studies (Kurniawan and Khardi, 2011; Zhou et al., 2019; Cui et al., 2022) focused on aircraft emissions during the LTO cycle. Unlike road, rail, and sea transportation, the flight process in the LTO cycle has prominent four-dimensional (4D) characteristics. For example, aircraft emissions have typical hourly temporal variations due to the impact of human activities. Moreover, the aircraft's unique three-dimensional (3D) flight trajectory (Koudis et al., 2017) makes it a distinctive 3D linear emission source. As a result, comprehensive spatial and temporal consideration is crucial for accurately calculating the pollutant emissions of aircraft in the LTO cycle.

For calculating pollutant emissions of aircraft in the LTO cycle, most of the current research is based on the ICAO standard method (Kurniawan and Khardi, 2011; Cui et al., 2022). ICAO stipulates that the LTO cycle is divided into four modes: take off, climb, approach, and taxi, reflecting that the standard operation time of each mode is 0.7, 2.2, 4, and 26 min, respectively (ICAO, 2011). However, unchanged running time is inconsistent with the actual aircraft operation process (Xu et al., 2020) because the running time of different modes in the LTO cycle is influenced by runway congestion (Badrinath et al., 2020) and MLH variations (Peace et al., 2006; Nahlik et al., 2016). Therefore, relying on the ICAO method may lead to high uncertainties. An alternative approach is to use accurate flight data, such as ADS-B data (Klenner et al., 2022; Zhang et al., 2022), which can significantly improve the accuracy of pollutant emission calculations. However, this method still has problems, such as difficulty obtaining actual aircraft data and limited application range. Therefore, multi-year hourly aircraft emission datasets that accurately reflect reality are still lacking. In air quality simulation, addressing the issue of pollutant emission inventory in the LTO cycle of aircraft in the spatial dimension is a significant challenge. Previous studies have primarily focused on the environmental impact of pollutant emissions from aircraft during the LTO cycle (Yim et al., 2015; Yang et al., 2018; Bo et al., 2019). However, most of these studies have allocated these emissions to the grid where the airports are without considering the altitude, longitude, and latitude of the emission locations. While this allocation method is suitable in rough-grid settings, using a finer grid to reflect aircraft emissions' environmental impact more accurately leads to more significant errors (Kumar et al., 1994; Arunachalam et al., 2011; Woody and Arunachalam, 2013). Therefore, considering the actual flight characteristics of aircraft is vital to obtaining more realistic spatial characteristics of aircraft pollutant emissions and improve the accuracy of air quality simulation. The impact of aircraft emission heights and horizontal position distribution modes on air quality varies widely, as demonstrated by various studies (Unal et al., 2005; Wolfe et al., 2016; Woody et al., 2016; Lawal et al., 2022). Zhang et al. (2023) conducted air quality simulations based on actual flight trajectories in the ADS-B data for typical regions. However, this method is limited by the availability of flights with ADS-B data and cannot be widely applied (Quadros, 2022). Consequently, there is still a lack of aircraft 4D emission inventory datasets in the LTO cycle that accurately reflect actual 3D flight trajectories and their dynamic nature over time.

As the world's second-largest aviation market (CAAC, 2021b), China contributes 13 % of global flight operations (Graver et al., 2020) and accounts for 7.8 % to 23.5 % of global aviation-related pollutant and carbon emissions (Ma et al., 2024; Teoh et al., 2024). Improving the accuracy of aviation emission estimates and enhancing temporal–spatial resolution in China can not only promote the green development of the Chinese aviation industry but also exert a far-reaching impact of global aircraft pollution mitigation. The period 2019–2023 is a unique period of the COVID-19 outbreak. Therefore, we have developed a 4D aircraft emission inventory (4D-LTO emission inventory dataset) for mainland China's takeoff and landing (LTO) cycle from 2019 to 2023. This inventory provides detailed and accurate emission calculations and flight trajectory recognition. It offers high spatial and temporal resolution, with a horizontal resolution of 0.03° × 0.03° and 34 layers of height resolution from 0 to 15 668 m. Our dataset and methodology are valuable resources for studying the temporal and spatial variations, air quality, and health impacts of aircraft emissions during the LTO cycle.

2 Methodology and data

Figure 1 illustrates the process of establishing the 4D-LTO emission inventory dataset, including the methods, the primary dataset, and the final output information. We developed the 4D-LTO emission inventory dataset in four steps:

1. accurately estimating the pollutant emissions of aircraft in the LTO cycle,

2. identifying the 3D flight trajectory of aircraft in the LTO cycle,

3. applying the spatial and temporal allocation method of the 4D-LTO emission inventory dataset,

4. comparing with the conventional spatial allocation method.

Figure 1 Schematic of the method used to develop a high-spatial- and high-temporal-resolution aircraft 4D emission inventory in the LTO cycle. It includes a detailed pollutant emission calculation method, flight trajectory recognition method (the horizontal recognition method takes the departure process as an example), and LTO cycle emission inventory allocation method. Publisher's remark: please note that the above figure contains disputed territories.

The 4D-LTO emission inventory dataset is a grid emission inventory dataset established by combining the aircraft emission calculation method and the flight trajectory identification method for the LTO cycle (described in Sect. 2.1 and 2.2.). The spatial–temporal allocation method is introduced in Sect. 2.3, and the comparison method is described in Sect. 2.4.

2.1 Aircraft LTO cycle emission calculation

The emission index, fuel flow, and running time of different flight modes are used to estimate the civil aircraft emissions in China based on the ICAO method (Kurniawan and Khardi, 2011; Bao et al., 2023). The calculation method is as in Eq. (1):

$$\begin{array}{}\text{(1)}& E_j=\sum _l\sum _m{\mathrm{EI}}_{l,m,j}\times {\mathrm{FF}}_{l,m}\times T_{l,m},\end{array}$$

where E_j is the emission (g) of pollutant j (including NO_x, HC, SO₂, CO, and PM) and EI_l,m is the emission index (g kg⁻¹) in m mode (take off, climb, approach, and taxi) of LTO l. ${\mathrm{FF}}_{l,m,j}$ is the fuel flow (kg s⁻¹) of pollutant j in m mode of LTO l, and T_l,m is the running time in m mode of LTO l.

The actual parameters of each flight should be used to calculate the emissions. However, complete data cannot be obtained due to problems such as incomplete data recording and recording errors. As a result, we have used different methods to approach the emission index, fuel flow (Sect. 2.1.1), and running time (Sect. 2.1.2 and 2.1.3) in Eq. (1) in the actual situation and estimate more accurate emissions.

2.1.1 Aircraft–engine matching

The aircraft's emission factor and fuel flow depend on its engine type, and the same aircraft type can be equipped with different types of engines. Thus, we collected as much detailed engine configuration information as possible for various aircraft types to improve the accuracy of the calculations. The matching method is divided into three steps:

1. We counted all aircraft types departing from or arriving in China from 2019 to 2023 using the flight information dataset from VariFlight and querying the aircraft type corresponding to each aircraft code.

2. We carefully counted China's airlines, civil aviation fleet in service information, point-type statistical engine number, type, and proportion through a comprehensive search in flight-associated dynamic query (VariFlight, 2024) and Civil Aviation Leisure Station (CALS, 2024).

3. We weighed the EI and FF of each aircraft type to obtain the value (Yang et al., 2018) using the information of all aircraft types and the proportion information of different engine types for each aircraft type combined with the emission index (EI) and fuel flow (FF) data of each engine type given in the ICAO Aircraft Engine Emissions Databank (EEDB, 2024).

The EI of an aircraft type in different modes was calculated as Eq. (2):

$$\begin{array}{}\text{(2)}& {\mathrm{EI}}_{i,m,j}=n_i\times \sum _k{\mathrm{EI}}_{k,m,j}\times P_{i,k},\end{array}$$

where ${\mathrm{EI}}_{i,m,j}$ is the emission index of aircraft type i in mode m (g kg⁻¹) of pollutant j (NO_x, HC, and CO), n_i is the number of engines fitted to aircraft type i, ${\mathrm{EI}}_{k,m,j}$ is the emission index of engine k in mode m of pollutant j (g kg⁻¹), and P_i,k is the proportion of aircraft type i equipped with engine k.

The FF of an aircraft type in different modes was estimated as Eq. (3).

$$\begin{array}{}\text{(3)}& {\mathrm{FF}}_{i,m}=n_i\times \sum _k{\mathrm{FF}}_{k,m}\times P_{i,k},\end{array}$$

where FF_i,m is the fuel flow of aircraft type i in mode m (kg s⁻¹), FF_k,m is the fuel flow of engine k in mode m (kg s⁻¹), and definitions of other parameters are similar to those used in Eq. (3).

In addition, the first-order approximation 3.0 (FOA3.0) (Wayson et al., 2009) method was used to recalculate the EI of PM, which is not included in EEDB. The emission factor of SO₂ is related to the sulfur content of jet fuel, so we used 3.868 g kg⁻¹ as the emission factor of SO₂ (GB6537, 2018). In summary, the references for the EI for different pollutants are shown in Table S1 in the Supplement.

2.1.2 Climb and approach time calculation

The daily maximum mixing layer height (MLH) serves as a key parameter for determining climb and approach modes of flight operations and varies with region and time. Given data accessibility constraints, we substituted daily maximum MLH with the daily maximum planetary boundary layer height (PBLH), which shares analogous dynamic characteristics. The three steps for calculating climb and approach times are as follows.

Different airport daily maximum PBLHs in 2019–2023 were obtained based on the Weather Research and Forecasting (WRF) model. The model parameter settings were described in our previous study (Wen et al., 2023).

The relationship between flight time and height was established. In our previous study (Zhou et al., 2019), the relationship between different airports in different months under the approach and climb mode was built based on Aircraft Meteorological Data Relay (AMDAR) data. AMDAR includes the aircraft's position (longitude, latitude, and altitude), speed, and associated meteorological parameters, which were collected by the aircraft navigation system. The recording intervals are set at 6 s for the first 60 s of the climb phase followed by once every 35 s thereafter and once every 60 s during the descent phase. The form of the relationship for climb and approach mode can be found in Sect. S1 in the Supplement. The R² (p<0.001) of the functional relationships of the climb and approach mode was above 0.93.

Each flight's actual climb and approach times from 2019 to 2023 were calculated based on the relationship between climb and approach mode mentioned above and the daily maximum PBLH at different airport.

2.1.3 Taxi-in and taxi-out time calculation

ICAO specifies the taxi mode's running time (taxi out 19 min; taxi in 7 min). However, the actual taxi time varies based on airport flight schedules during actual operation, and using a fixed time can lead to emission calculation uncertainty. Therefore, the actual taxi time data were used to calculate the aircraft's taxi emissions accurately. The actual taxi time data were obtained from VariFlight based on the information of the ADS-B system, which is recognized by researchers as a reliable data source (Klenner et al., 2022; Zhang et al., 2022; Teoh et al., 2024). Since not all aircraft record the actual taxi time and the actual taxi time is not publicly available (Table S2), this study collects all available taxi time data, with coverage rates ranging from 47.9 % to 67.0 % during 2019–2023, which are summarized in Table S2. They can represent the taxiing conditions of aircraft at different airports with different operating scales. The missing taxi time was supplemented based on the hourly airport difference relationship model between taxi time and aircraft number of schedules. The functional relationship between the number of aircraft on schedule and the taxi time is as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }{T}_{\mathrm{taxi}}=\mathrm{\Delta }T\times \mathrm{Ns}+T_{\mathrm{0}},\text{(5)}& {\displaystyle }\mathrm{\Delta }T=u\cdot e^{vN},\text{(6)}& {\displaystyle }{T}_{\mathrm{0}}=o\cdot e^{dN},\end{array}$$

where T_taxi is the taxi-out (in) time of each flight (s) and ΔT is an increase in taxi time per Ns (s per aircraft). T₀ is the initial taxi time (s), and Ns is the number of aircraft on schedule in an hour. N is the annual average aircraft departure/arrival number for each hour; u, v, o, and d are the airport-specific constants.

The hourly airport difference relationship model between T_taxi and Ns at different airports was used to update the database from 2019 to 2023. The performance of the taxi time calculating model for different airports is shown in Fig. S1 and Table S3, taking 2023 as an example. In addition, in this study, the Beijing Capital International Airport (PEK) is chosen as a case to test the performance of taxi time model (Fig. 2a and b) in diverse flight situations (e.g., high-density scenarios), due to the centralized terminal layout and relatively frequent ground congestion (Liu et al., 2024). Taking 12:00 UTC+8 from 2019 to 2023 for the Beijing Capital International Airport (PEK) as an example, Fig. 2a, b represent the comparative verification of function relationships for taxi-in and taxi-out modes in different years. We observed a strong correlation between taxi time and the number of scheduled aircraft, regardless of whether it is taxi in or taxi out. The significance level (p<0.001) indicates a strong relationship. The R² for taxi-out mode ranges from 0.87 to 0.98, and for the taxi-in mode, it ranges from 0.96 to 0.99. The model has a good effect on taxi-in or taxi-out mode at different years, indicating that the model reflects the real taxi time variation.

If the relationship between taxiing time and the number of aircraft scheduled cannot be fitted to a certain time due to lack of records, Table S4 presents the exponential relationship of ΔT and T₀ at different years. In addition, Fig. 2 also provides the 5-year exponential relationship of ΔT and T₀, which could be a reference for the other study with no fitting data. We also calculated the coefficient of variation (CV; 30.4 % for taxi-out and 10.4 % for taxi-in operations) between the ΔT and T₀ estimation result from the 5-year model and the specific-year model for a representative flight number of 20 (common across all study years). Compared with the actual taxi time, the estimation error in the model result (11.4 % for taxi-in mode and 20.4 % for taxi-out mode) is lower than the result based on the fixed ICAO standard taxi time (27.8 % for taxi-in mode and 22.0 % for taxi-out mode). ΔT and T₀ estimation models are only used in the situation when the ΔT and T₀ cannot be counted due to a lack of records.

Figure 2 The linear function relationship between taxi time (T) and the number of aircraft on the schedule at each hour (Ns), the exponential function relationship between ΔT and T₀, and the annual average departure aircraft number for each hour (N). (a) Taxi out. (b) Taxi in.

2.2 Aircraft emission 3D trajectory identification

This study is divided into two steps to identify the 3D spatial location of aircraft emissions in China during 2019–2023: (1) the flight altitude identification (Sect. 2.1.1) and (2) each flight's horizontal trajectory identification (Sect. 2.1.2).

2.2.1 Flight altitude identification

The relationship between flight height and time (Zhou et al., 2019), which is introduced in Sect. 2.1.2, was used to identify the altitude at different times for each LTO cycle. The daily maximum PBLH was used to identify the maximum height of each LTO cycle at different airports. When the taxi in and out are 0 m, the takeoff is from 0 to 152 m (ICAO). The climb is from 152 m to PBLH, and the approach is from PBLH to 0 m. By integrating the altitude information with the emission inventory data established in Sect. 2.1, we were able to further vertically stratify the pollutant emission inventory during the LTO cycle.

2.2.2 Flight horizontal trajectory identification

We established the flight trajectory database of each airport in China based on the density-based spatial clustering of applications with noise (DBSCAN) algorithm. Moreover, each flight's trajectory was identified based on the relative position of the departure airport, the arrival airport, and the wind direction. The clustering method is used to screen out flight trajectories with similar characteristics from a large amount of actual flight data. This approach helps determine the grid location of aircraft emissions (Gariel et al., 2011; Bombelli et al., 2017). DBSCAN is a density-based clustering algorithm widely used in machine learning and data mining (Chen et al., 2021b; Tekin and Sarı, 2024). For the transportation industry, it is used for the identification research of road traffic, ship, and aircraft trajectories (Gui et al., 2021; Deng et al., 2023; Li et al., 2023). The DBSCAN algorithm belongs to unsupervised learning, and the initial value setting does not significantly affect the clustering results (Ventorim et al., 2021). As a result, the DBSCAN algorithm is well suited for flight trajectory clustering processing with unclear information, such as the number of clusters and distribution characteristics (Murça et al., 2018; Giovanni et al., 2024).

Before clustering, flight trajectory data belonging to the LTO cycle should be extracted from a vast amount of information in AMDAR. First, the climb and approach modes in the LTO cycle are screened according to the ascending and descending symbols in AMDAR information. Second, each flight trajectory is divided based on airport ownership according to the airport's location. Finally, the horizontal position information (time, longitude, and latitude) of each flight trajectory in the climb and approach modes of different airports is obtained as the input information for flight trajectory clustering.

The DBSCAN algorithm relies on two input parameters, the minimum number of samples (MinPts) and the distance threshold (ε), to cluster the data space based on three basic concepts: directly density-reachable, density-reachable, and density-connected (Sander et al., 1998). MinPts determines the minimum number of points required to form a dense region, while ε specifies the maximum distance between two points to be considered to be within the same neighborhood.

DBSCAN is good at calculating the distance between points, but it is difficult for DBSCAN to process the flight trajectory with the time attribute in this study (Chen et al., 2021a). Therefore, we use the Euclidean norm to compute the distance between the two sets of flight trajectories. The premise of using the Euclidean norm is to keep the time intervals of each set of flight trajectories the same. However, the time interval of each flight trajectory sequence is not the same because of each flight's trajectory difference and recording delay. As a result, we conducted unified processing of each departure and arrival trajectory using the resampling method. Sampling points that are too low and too high make the location feature information unclear and increase the computational complexity of clustering processing, respectively. Based on all actual flight data from 2019–2023, during the LTO cycle, departure was within 480 s and arrival was within 1200 s. To comprehensively consider the recording intervals of AMDAR data, the uniformity across departure and arrival phases, and computational complexity, we set the sampling points of each trajectory to 25.

MinPts selection was performed as follows: since the average number of sample trajectories varies in different airports, the MinPts must be determined separately for various airports. For each airport, the MinPts values are, respectively, taken to be in the range of 6 to 10, and the clustering effect is observed, from which the appropriate MinPts values are selected. ε selection was performed as follows: the method k-distance (Garg et al., 2020) graph method was used to select the appropriate ε. The k-distance curve first calculates the distance between each trajectory in the data and the trajectory with the nearest k and then arranges the k-distances of all trajectories in descending order and draws the curve. Moreover, k values for different airports are the same as MinPts, and the ε values are based on the apparent inflection point in the k-distance curve.

Figure S2 shows the overall performance of the trajectory clustering model for different airports. In addition, in this study, the Shanghai Pudong International Airport (PVG) is chosen as a case study to test the performance of trajectory cluster (Fig. 3) under the normal flight trajectory scenario as well as the deviations in flight trajectories due to the crosswinds and typhoons, which is the challenge for the robustness of the trajectory cluster algorithm (Wang et al., 2017; Xu et al., 2020). This study evaluates the core flight trajectory of departure and arrival since the flight trajectory is a series of latitude and longitude information with time series characteristics. This approach uses the DBSCAN clustering method by splitting the flight trajectory of departure and arrival into two directions, latitude and longitude, considering the three indices of R and MAE. In longitude, the correlation between identified core trajectories and actual trajectories is more significant than 0.80 (0.865–0.992 for departure; 0.811–0.997 for arrival). In latitude, the correlation between identified core trajectories and actual trajectories is more significant than 0.94 (0.947–0.995 for departure; 0.941–0.992 for arrival). The identified core trajectory is consistent with the actual flight trajectory, indicating that the core trajectory can reflect the actual flight situation. In longitude, the MAE between identified core trajectories and actual trajectories is less than 0.05° (0.01–0.02° for departure; 0.02–0.05° for arrival). In latitude, the MAE between identified core trajectories and actual trajectories is less than 0.05° (0.01–0.02° for departure; 0.02–0.05° for arrival). Although the clustering results are uncertain, they can still provide vital information for the 3D grid location of aircraft emissions.

Figure 3 (a) Departure flight trajectories of different clusters. (b) Consistency in longitude of departure between the core and actual flight trajectories. (c) Consistency in latitude of departure between the core and actual flight trajectories. (d) Arrival flight trajectories of different clusters. (e) Consistency in the longitude of arrival between the core and actual flight trajectories. (f) Consistency in the latitude of arrival between the core and actual flight trajectories.

The airports with multiple runways assign a suitable runway for each flight based on the relative location of the departure and arrival airport (Yin et al., 2022; Sekine et al., 2023). This decision is made considering the need for aircraft to operate against the wind (CMA, 2015; CAACNEWS, 2019) as per the Chinese Meteorological Administration. Therefore, these flight characteristics were combined to identify the horizontal trajectory of each flight in the LTO cycle (horizontal position recognition in Fig. 1). First, all the flight trajectory clusters corresponding to the departure/arrival airport are selected from the flight trajectory database obtained by the DBSCAN method. Second, trajectories from the runway are chosen to be close to the target airport. Third, the aircraft takes off against the wind principle, selecting trajectories and the side of the runway based on the wind direction information at the moment of the departure from/arrival at the airport. Finally, the final trajectory is selected by the target direction being the opposite of or the same as the takeoff direction.

2.3 Temporal and spatial identification of 4D emission inventory

Gridded emission information is often required for air quality and climate simulation models or refined prevention and control of pollutants. Therefore, the obtained aircraft pollutant emission inventory of the LTO cycle in China during 2019–2023 was processed into hourly 3D-grid pollutant emission data with a horizontal resolution of 0.03° × 0.03° and a height resolution of 34 layers from 0 to 15 668 m.

2.3.1 Aircraft emission temporal allocation

The 4D-LTO emission inventory dataset has an hourly temporal resolution. According to Eq. (1), the emissions of each pollutant in different modes of each LTO cycle are calculated separately, and the emission for each hour of the LTO cycle at different airports is the sum of the pollutant emissions generated by all departure and arrival at that airport during that hour. In addition, the daily, monthly, and yearly total emissions are the sum of all LTO cycles of that day, month, and year to further analyze the temporal variation in pollutant emission.

2.3.2 Aircraft emission spatial allocation

For the horizontal resolution, most airport runways are approximately 3–4 km (CAAC, 2021a) in length, and certain pollutants (such as CO) are predominantly emitted during taxiing, i.e., on the runway. Overall, 0.03° × 0.03° is capable of reflecting the horizontal distribution characteristics of aircraft emissions. In addition, 0.03° × 0.03° is also a common resolution for air quality models. Therefore, the horizontal resolution of the 4D-LTO emission inventory is 0.03° × 0.03°, with a latitude and longitude range of 3.40–53.56° N and 73.44–135.09° E, respectively.

For the altitude resolution, while ICAO defines the LTO cycle with a fixed mixing layer height (915 m), in reality, the mixing layer height varies significantly with region and time, leading to variations in the altitude range of the LTO cycle. Therefore, to better reflect the vertical distribution of aircraft emissions above 915 m during the LTO cycle, this study sets the altitude range to from 0 to 15 668 m. In addition, to ensure that the emission inventory can be effectively used in air quality models, this study uses the air quality model commonly used 35-layer sigma stratification strategy (Wolfe et al., 2016). Therefore, the altitude resolution was divided into 34 layers from 0 to 15 668 m (0.0–38.3, 38.3–76.7, 76.7–115.3, 115.3–154, 154–231.8, 231.8–310.3, 310.3–389.3, 389.3–469, 469–549.3, 549.3–630.3, 630.3–711.9, 711.9–794.2, 794.2–960.7, 960.7–1130.1, 1130.1–1302.3, 1302.3–1477.6, 1477.6–1656.0, 1656.0–1929.7, 1929.7–2211.1, 2211.1–2599.3, 2599.3–3107.2, 3107.2–3643.1, 3643.1–4210.5, 4210.5–4813.9, 4813.9–5458.5, 5458.5–6151.2, 6151.2–6900.4, 6900.4–7717.4, 7717.4–8617.3, 8617.3–9621.2, 9621.2–10 759.7, 10 759.7–12 080.6, 12 080.6–13 664.8, and 13 664.8–15 668 m).

The 4D-LTO emission inventory dataset was processed by first identifying the emission information of each flight into a 3D grid using latitude, longitude, and altitude information. Then, the emissions of all flights within the same hour were summarized.

2.4 Comparison of our dataset with the previous dataset

Our dataset was compared with the spatial allocation methods commonly used in previous studies. (1) Other studies typically assign aircraft emissions in the LTO cycle according to the standard altitude for each mode as defined by ICAO (Mokalled et al., 2018; Bo et al., 2019; Wang et al., 2023; Zhang et al., 2023). (2) The conventional horizontal distribution method for aircraft emissions in the LTO cycle assumes that aircraft emissions are radially distributed (Lawal et al., 2022). The Federal Aviation Administration (FAA, 2024) recommended the standard climb rate of 200 ft per nautical mile. Therefore, the standard climb rate and ICAO standard altitude determine the horizontal distribution of aircraft emissions around the airport. The running time, altitude, and horizontal range of each mode defined by ICAO are shown in Table 1.

Table 1 The running time and altitude range of each mode defined by ICAO.

2.5 Uncertainty calculation

The uncertainty of the 4D-LTO emission inventory dataset is mainly divided into emission calculation uncertainty and spatial location identification uncertainty. This study assumes that the uncertainty in all input parameters follows a normal distribution.

When calculating the emission uncertainty, this study comprehensively considers the uncertainty in EI, FF, and T. The EI and FF are weighted based on the engine data from the EEDB and the engine proportion data for different aircraft types. Therefore, the standard deviation of EI or FF was calculated using Eq. (7):

$$\begin{array}{}\text{(7)}& \mathit{\sigma }=\sqrt{\sum _k(x_k-\overline{x}{)}^{\mathrm{2}}\times P_k},\end{array}$$

where σ represents the standard deviation of EI or FF, k represents the engine type, x_k represents the EI or FF value of engine k, $\overline{x}$ represents the weighted average of EF or FF, and P_k represents the proportion of engine k.

The climb and approach time is obtained using the relationship between flight time and flight height (Zhou et al.,2019). Therefore, the standard deviation of the climb and approach time is the combination of the standard deviation of function fitting parameters a, b, and c. The taxi-in/taxi-out time is calculated using Eq. (4). Therefore, the standard deviation of the taxi-in and taxi-out time is the combination of the standard deviation of a function fitting parameter ΔT and T₀. This study uses the Monte Carlo sampling method to obtain the 95 % prediction interval of the emission for different pollutants with 20 000 samples.

The spatial uncertainty during the LTO cycle includes the uncertainty in horizontal and altitude positions. The standard deviation of the horizontal position is calculated by the error distribution between the flight trajectory clustering result and the actual flight trajectory. The standard deviation of the altitude position is the combination of the standard deviation of function fitting parameters a, b, and c. This study employs the Monte Carlo method to quantitatively assess the uncertainty of spatial location identification for each hour, with uncertainty ranges derived from 20 000 Monte Carlo simulations at a 95 % prediction interval.

3 Results and discussion

3.1 Total aircraft emissions in the LTO cycle

In 2023, the total emissions of five types of pollutants in the LTO cycle of aircraft in China are as follows: HC is 3.2 Gg, CO is 46.1 Gg, NO_x is 62.3 Gg, PM is 1.1 Gg, and SO₂ is 18.4 Gg, as shown in Fig. 4a. The annual emission of different pollutants in 2023 was 82.9 % (HC)–94.1 % (NO_x) in 2019. However, before 2022 (the last year impacted by COVID-19 and the most affected year), emissions of various pollutants averaged 34.7 %–42.8 % of 2019. At the end of COVID-19, the 2023 recovery in aircraft emissions shows that the pandemic did not have an irreversible impact on aircraft activities and that emissions from aircraft activity will continue to grow (Teoh et al., 2024). Emissions of pollutants from aircraft, such as NO_x and PM_2.5, are known to cause respiratory and cardiovascular issues (Boningari and Smirniotis, 2016; Hu et al., 2022; Hou et al., 2024). Therefore, it is essential to pay attention to the growing trend of aircraft activities in order to anticipate and address their potential health impacts.

Figure 4a shows the main emission contribution of HC and CO came from taxi mode (94.6 % for HC; 91.5 % for CO) because HC and CO are mainly produced by incomplete fuel combustion, taking 2023 as an example. A large amount of HC and CO is created because the engine's thrust in taxi mode is minimal and the operation time is long (EPA, 1981). The climb is the main NO_x emission stage (42.1 %). The takeoff with the shortest running time contributes to the second-largest NO_x emission (25.7 %). The taxi with the longest running time contributes the most minor NO_x emission of 12.4 %, indicating that the emission factor of NO_x is highly correlated with the aircraft engine's thrust (Stettler et al., 2011). Although the engine runs for a long time, the NO_x emission during taxi mode with a slight thrust is still lower than during the takeoff stage, with the engine running for a short time but at nearly full thrust. For PM and SO₂, the emission contribution ratio is similar to the running time of each mode, and the taxi mode with the longest running time contributes 33.1 % of PM and 35.1 % of SO₂. The climb mode contributes 28.4 % of PM and 26.6 % of SO₂, the approach mode contributes 25.7 % of PM and 26.4 % of SO₂, and the takeoff mode with the shortest running time contributes 12.7 % of PM and 11.9 % of SO₂. From 2019 to 2023, among various aircraft types, B738, A320, and A321 were the top three pollutant emissions (Fig. 4b). The top three aircraft types contributed 64.1 % of NO_x emissions in 2019, taking NO_x emissions as an example. However, during the COVID-19 period (2020–2022), the contribution of the top three aircraft types reached 70.3 %–70.9 %. At the end of the pandemic impact in 2023, the contribution of the top three aircraft types reversed to the 2019 level (55.0 %). During the COVID-19 pandemic, many aircraft types ceased operation, including F50, E145, other regional aircraft, A306, A340, and other wide-body aircraft types, increasing the proportion of the first three types. As the impact of COVID-19 gradually diminished, the discontinued models resumed operation, and the emission proportion of the first three models returned to normal.

Figure 4 (a) Total aircraft pollutant emissions of the LTO cycle from 2019 to 2023. (b) Proportion of NO_x emissions in different aircraft types from 2019 to 2023.

3.2 Temporal variation in aircraft emissions in the LTO cycle

Figure 5a shows the changes in aircraft emissions during the LTO cycle from 2019 to 2023, encompassing the period before, during, and after the COVID-19 pandemic. The baseline year for analysis is 2019, unaffected by COVID-19, and represents regular aircraft activity.

As can be seen from Fig. 5a and Table S5, from 20 January to 13 February 2020, aircraft activity rapidly dropped to the lowest point owing to the impact of COVID-19, showing that the number of LTO on 13 February 2020, was 84.8 % lower than the same period in 2019. In the following months, aircraft activity slowly recovered, returning to the 19-year level in October. As the COVID-19 situation in China entered a recurrent period, from 2021 to the beginning of 2023, the activity of aircraft fluctuated, reflecting five low points (12 February 2021, 12 August 2021, 9 November 2021, 4 April 2022, and 19 November 2022). As the effects of COVID-19 faded from early 2023 on, aircraft activity gradually returned to 2019 levels. When the impact of COVID-19 ended, abnormal growth is noted in aircraft activity. In May 2021, the number of aircraft LTO increased rapidly compared to the same period in 2019. However, during the same period in 2019, aircraft activity showed a downward trend. From July to October 2023, the number of LTOs exceeded the same period in 2019. This phenomenon occurs because people with unfulfilled travel needs are inclined to engage in revenge tourism following prolonged COVID-19 lockdowns, resulting in increased aircraft activity and a sudden increase in emissions in the short term.

Figure 5 (a) Daily variation in pollutant emissions and the number of LTO cycles from 2019 to 2023. (b) Annual hourly variation in pollutant emissions and the number of LTO cycles from 2019 to 2023.

Based on Fig. 5b, the emission of various pollutants in 2019–2023 varies slightly in hours, with higher daytime (06:00–23:00 UTC+8) and low nighttime (00:00–05:00 UTC+8) values, with the minimum at 04:00 UTC+8 and the maximum at 13:00 UTC+8. The most significant difference in the number of LTOs and the emission of pollutants between each hour over the 5 years occurred at 04:00 and 13:00 UTC+8 in 2022 (251 % of LTO, 230 % of HC, 244 % of CO, 229 % of NO_x, 249 % of PM, and 234 % of SO₂). The difference in pollutant emissions between 2019 and 2023 at each hour shows recovery in 2023. The proportion of LTO numbers that recovered to 2019 levels in the nighttime was 34.1 % lower than during the daytime. Additionally, the recovery rate of the five pollutants' emissions in the nighttime was 39.1 %–44.4 % lower than in the daytime, indicating that resumed aircraft activity was significantly better during the day than at night.

3.3 4D characteristics of aircraft emissions in the LTO cycle

During the LTO cycle, HC and CO were predominantly emitted during taxi mode (Yang et al., 2018). Consequently, HC and CO emissions are distributed in the first layer of the grid where the runway is located. NO_x is an important contributor to overall aircraft emissions and has a significant impact on air quality (Zhang et al., 2023). Furthermore, the spatial distribution of PM and SO₂ emissions from aircraft is similar to that of NO_x. In summary, this study mainly analyzes the spatial distribution of NO_x emissions.

This study calculates hourly aircraft emissions in LTO cycles at various airports in China during 2019–2023 based on the combining emission calculation method and flight trajectory recognition method, establishing a 4D aircraft NO_x emission inventory (hourly, 0.03° × 0.03° × 34) of the LTO cycle in China (Figs. 6 and 7).

Figure 6a signifies the horizontal distribution of yearly NO_x emissions in prefecture-level cities and airports during 2019–2023. Compared with 2019, emissions in most regions affected by COVID-19 decreased significantly during 2020–2022. Notably, aircraft emissions of prefecture-level cities experienced an average reduction of 43.1 % in 2022. As the COVID-19 impact ended in 2023, aircraft emissions of prefecture-level cities recovered, with an average increase of 5.07 %. Although the aircraft emissions of the LTO cycle in prefecture-level cities fluctuated from 2019 to 2023, airport emissions in Beijing, Shanghai, Guangzhou, and Chengdu were the top four, accounting for 24.4 % (2022)–32.2 % (2023) of national emissions.

Between 2019 and 2023, the number of airports in China increased from 237 to 257 at an average annual growth of 5. However, the newly operated airports significantly increased aircraft emissions in a prefecture-level cities. For example, due to the operation of the Chengdu Tianfu International Airport (TFU), aircraft NO_x emissions in Chengdu (4.2 Gg) were 32.3 % in 2023, which is higher than in 2019. In addition, Chengdu's aircraft NO_x emissions were 8.4 %–14.3 % higher than Guangzhou's during 2021–2023, while in 2019, Chengdu's NO_x emissions were 21.5 % lower than Guangzhou's when the TFU airport did not start operations. The newly operated airports can also affect the original airport in a prefecture-level cities. Taking airports in Beijing as an example, PEK airport's annual aircraft NO_x emissions (8.1 Gg) were 101 %, ranked second in 2019, which is higher than those of CAN airport (4.0 Gg). However, with PKX airport's operation, PEK airport emissions significantly decreased. In 2023, PEK airport's emissions recovered to 54.7 % of the original, while the total emissions of Beijing recovered to 81.4 %. In addition, emissions from the PEK airport in 2023 were only 15.7 % higher than those from CAN, indicating that the newly operated PKX airport has reduced the emission pressure on PEK airport.

Taking airports in Beijing and surrounding areas in January 2023 as an example, Fig. 6b demonstrates the grid horizontal distribution of aircraft NO_x emissions in the LTO cycle. The horizontal distribution characteristics of aircraft emissions in the LTO cycle are influenced by the distance along the runway and how they spread, indicating that emissions are concentrated in the direction of the runway near the airport. With the increase in flight distance, the emissions caused by aircraft are dispersed. Aircraft emissions during the LTO cycle are widely distributed around the airport and not even represented by a rough grid (e.g., 0.36° × 0.36°). The elongated distribution characteristics of aircraft emissions indicate that evenly allocating emissions around the airport will cause significant uncertainty. Figure 6c shows the differences in aircraft emissions at various airports and times between 00:00 and 20:00 UTC+8 on 3 January 2023 at a 4 h interval. This phenomenon indicates that the horizontal distribution characteristics of aircraft emissions vary significantly at different hours and airports. As a result, the refined aircraft emission inventory on the LTO cycle conforms to the time-by-hour spatial distribution characteristics of aircraft, better reflecting the actual situation of aircraft emissions, which is of great significance for accurately assessing aircraft environmental impact in the LTO cycle.

Figure 6 (a) Horizontal distribution of yearly NO_x emissions in prefecture-level cities and airports during 2019–2023. (b) Horizontal distribution of NO_x emissions at airports in Beijing and surrounding areas in January 2023. (c) Horizontal distribution of NO_x emissions at airports in Beijing and surrounding areas for different hours in January 2023. Publisher's remark: please note that the above figure contains disputed territories.

Figure 7 (a) Height distribution characteristics of NO_x emissions in the LTO cycle. (b) The maximum altitude layer for different months and the corresponding altitude range. (c) The emissions above and below the height of 915 m in different months and the proportion above 915 m. (d) Height distribution characteristics of NO_x emissions from the LTO cycle at different months.

Figure 7a uses the annual NO_x emissions in 2023 to demonstrate the height distribution of aircraft emissions in the LTO cycle. In general, the NO_x emission of aircraft in the LTO cycle decreases with the increase in altitude. Moreover, the emission per unit altitude significantly decreases between layers 1 and 2 and between layers 4 and 5 due to the different flight altitude ranges in various modes in the LTO cycle. Emissions from layer 1 (0–38 m) include the entire taxi mode as the takeoff mode and approach mode, with the maximum unit height NO_x emissions (0.32 Gg m⁻¹). Emissions from layers 2 to 4 (38–154 m) include a part of takeoff mode and approach mode, with the unit height NO_x emissions of 0.11–0.12 Gg. From layer 5, each layer's NO_x emissions (≤ 0.04 Gg m⁻¹) include the part of the climb and approach modes. As the emission height increases, the emissions of NO_x gradually decrease. The reduction rate gradually increases before layer 14 and decreases after layer 14, indicating that the unit height emissions of each layer above the 14th layer have little difference. In addition, there are significant differences in the height distribution characteristics of emissions in the LTO cycle at different months. Figure 7b shows that the maximum emission height in the LTO cycle can reach the 23rd layer (3107–3643 m, November and December) and the 26th layer (5459–6151 m, June) of 34 layers (0–15 668 m). The maximum aircraft emission height in the LTO cycle can reach 4544 m above the ICAO-defined maximum altitude of 915 m due to MLH variation across 12 altitude levels. Figure 7c illustrates that the NO_x emissions above the ICAO standard height (∼ 915 m) in different months accounts for an average of 24.6 % (9.9 %–37.5 %) in the LTO cycle. This result indicates that the ICAO method does not account for a significant portion of emissions during the entire LTO cycle. Based on previous studies (Köhler et al., 2008; Lee et al., 2013; Yim et al., 2015; Zhang et al., 2023), high-altitude emissions can significantly impact ground-level air quality through atmospheric transport and chemical reactions. When assessing emissions during the LTO cycle and their impact on air quality and health, we must fully consider the contribution of emissions above 915 m. Therefore, using the ICAO fixed flight height introduces considerable uncertainty when calculating the aircraft emission during the LTO cycle and assessing its environmental impact.

3.4 Comparison with the previous allocation method

Figure 8 uses the NO_x emissions in January 2023 to show the differences between the 4D-LTO emission inventory and the LTO emission divided in previous studies (Mokalled et al., 2018; Bo et al., 2019; Lawal et al., 2022; Wang et al., 2023; Zhang et al., 2023) in terms of height distribution (Fig. 8a–b) and horizontal distribution (Fig. 8c–e). Figure 8a and b represent noticeable differences in emissions at different layer heights. Two statistical measures, mean absolute error (MAE) and mean absolute percentage error (MAPE), were employed to quantify these differences. (1) The first of two components of the allocation error of the ICAO method is in the range of 154–961 m; the ICAO method overestimates the emissions by 63.4 %, and the emission difference between different layers is 65.7–219.3 Mg. The difference increases with a rise in height. (2) The second component lies within the range of 961–4211 m, the ICAO method missed 283.0 Mg of emissions, and the difference decreases with an increase in height (0.0–125.7 Mg). Figure 8b uses PEK, PVG, and CAN to demonstrate the emission height changes between different airports. Different airports' overestimation and missing zones are similar to the height distribution of total NO_x emissions. However, the ICAO method misses emissions above 961 m differently for different airports (1.9 Mg for CAN, 28.0 Mg for PEK, 5.4 Mg for PVG), and the ICAO method overestimates emissions within 154–961 m differently at various airports (3.4–11.8 Mg for CAN, 3.4–10.2 Mg for PEK, 1.5–10.1 Mg for PVG). Compared with the dataset based on the ICAO method, our 4D-LTO emission inventory dataset can more accurately represent the height distribution of actual aircraft emissions.

In the example of airports in Beijing and surrounding areas, Fig. 8c and d demonstrate that our 4D-LTO emission inventory dataset outperforms the dataset based on the previous radial allocation method, showing an apparent misallocation of emissions. Figure 8e quantifies the differences in the horizontal distribution between two emission inventory datasets. Based on the previous radial allocation method, the dataset misallocated 242.7 Mg of emissions in the misallocation zone. Among them, 17.2 Mg of emissions were missing (3.0 Mg for PEK, 13.5 Mg for PKX, 0.2 Mg for TSN, 0.5 Mg for TVS), and 225.5 Mg of emissions were assigned to the wrong grid (122.8 Mg for PEK, 73.7 Mg for PKX, 25.8 Mg for TSN, 3.2 Mg for TVS). In the non-misallocation zone, the dataset based on the previous radial allocation method underestimates 41.9 % of emissions (46.5 % for PEK, 37.8 % for PKX, 32.9 % for TSN, 60.6 % for TVS). Compared with the dataset based on the previous radial allocation method, our 4D-LTO emission inventory dataset can better reflect the horizontal distribution of actual aircraft emissions.

Figure 8 Comparison of horizontal and height distributions of NO_x emissions in January 2023, (a) NO_x emission differences at different heights between two datasets (Mg), (b) NO_x emission distribution at different layers between two datasets (Mg), (c) distribution of NO_x emissions based on previous radial allocation method, (d) distribution of NO_x emissions in the 4D-LTO emission inventory dataset, and (e) NO_x emission differences at different horizontal grids between two emission inventory datasets.

3.5 Uncertainty analysis

Taking the year 2023 as an example, this study estimated the hourly uncertainty ranges for different emissions in various airports throughout the year, with the average of the uncertainty ranges for NO_x, CO, HC, PM, and SO₂ being [−17 %, 17 %], [−16 %, 16 %], [−6 %, 6 %], [−8 %, 8 %], and [−8 %, 8 %], respectively, and their standard deviations being ±11 %, ±10 %, ±3 %, ±4 %, and ±4 %, respectively.

This study also estimated the hourly average 95 % prediction intervals for latitude and longitude of departure and arrival flights at various airports, with the average of [−0.02°, 0.02°] (longitude) and [−0.02°, 0.02°] (latitude) for departure and [−0.06°, 0.06°] (longitude) and [−0.05°, 0.05°] (latitude) for arrival, respectively, and their standard deviations being ±0.01°, ±0.01°, ±0.01°, and ±0.01°, respectively. In addition, the result showed that the hourly 95 % prediction intervals of climb and approach altitude at different airports are within an average of [−78 m, 78 m] and [−134 m, 134 m], respectively, and that their standard deviations are ±49 m and ±95 m, respectively.

3.6 Advantage and limitation

The previous studies (e.g., Zhang et al., 2022, and Teoh et al., 2024) contributed a lot to the improvement of the emission estimation using ADS-B data, laying a solid foundation for further assessing the impact of aircraft emissions on air quality. In percentages, the ADS-B data do not fully cover all flights. It is hard to identify the flight trajectories of those flights not covered by ADS-B. Furthermore, using a fixed height range (915 m) for the LTO cycle introduces errors in the calculation of pollutant emissions during the LTO cycle. The other emission inventories (e.g., EDGAR, EMEP, and AERO2k), are also important sources of aviation emission data. These emission inventories mainly rely on the flight schedule information, which mainly calculate the total emissions of countries, and does not reflect the detailed four-dimensional emission characteristics.

Based on the emission index and activity level data of each mode during the LTO cycle, we calculate the emissions of each flight by the bottom-up method and give the hourly and 3D spatial distribution of aircraft emissions. It is important information for further assessing aircraft emissions during the LTO cycle and their impact on air quality.

Our 4D-LTO emission inventory dataset reflects the actual spatial and temporal and can be used to accurately assess the air quality impact of aircraft in the LTO cycle, but it has several limitations due to data and technical restrictions. (1) According to our investigation (Airbus, 2025; Aircraft Commerce, 2025), most aircraft types do not update engine. Some aircraft types (e.g., the aircraft A is allowed to be equipped with engines A and B) may change engine configuration proportion during 2019–2023. Given the unavailability of the annual variation of engine configuration for each aircraft type in the existing datasets, this study used the latest proportion data of each aircraft type in different years. (2) The certified engine emission indices derived from the engine manufacturers and reported in the ICAO failed to consider the life expectancy of an aircraft and meteorological conditions. This may result in errors between the fuel consumption and emissions estimated using these recommend parameters and real-world conditions. Therefore, future research should be conducted on the dynamic emission factors based on the machine age and flight conditions. (3) Our dataset was obtained from a subset of available flight data and near-real flight data based on a model built using real data. This may result in errors between our dataset and real-world conditions. However, this issue will be addressed as real-world data become more widely available.

4 Data availability

The 4D-LTO emission inventory dataset in China from 2019 to 2023 presented in this study is freely available at https://doi.org/10.5281/zenodo.13908440 (Lang et al., 2024).

5 Conclusions and implication

This study establishes China's 4D-LTO aircraft emission inventory dataset during 2019–2023 by combining accurate and generalizable emission methods and flight trajectory identification methods. The actual taxi time is used, and the supplementary value is obtained through the 5-year validation (for PEK, R²=0.87–0.99) airport hourly difference relationship between the taxi-in/taxi-out time and the number of aircraft scheduled. Moreover, the climb and approach time and the attitude of each flight are updated using the MLH and the airport monthly difference relationship between the flight altitude and time of the climb/approach mode. Finally, the DBSCAN clustering method (for PVG, R=0.865–0.995 and MAE = 0.01–0.02° in departure and R=0.811–0.997 and MAE = 0.02–0.05° in arrival) is used to obtain the flight trajectory database of each airport based on the massive number of actual flight trajectory data. Then, the flight trajectory of each flight is identified by the wind direction and the relative position of the departure and arrival airport. The data show that the impact of COVID-19 reduced the LTO number to 73.1 % in 2020, 77.6 % in 2021, and 48.7 % in 2022, compared to 2019. However, in 2023, the emissions of different pollutants quickly bounced back to 82.9 %–94.1 % of the 2019 levels, resulting in HC, CO, NO_x, PM, and SO₂ emissions of 3.2, 46.1, 62.3, 1.1, and 18.4 Gg, respectively.

Taxi is the most crucial emission stage of HC and CO (94.6 % and 91.5 % of the emission of the entire LTO cycle), and climb is the primary emission stage of NO_x (42.1 %). We also find that takeoff with the smallest operation time contributes the second-largest emission of NO_x (25.7 %). Moreover, B738, A320, and A321 are the top three aircraft types that emit pollutants. During the COVID-19 period (2020–2022), the contribution of the top three aircraft types reached more than 70 %.

Due to the impact of COVID-19, aircraft emissions in the LTO cycle fluctuate from 2019–2023. After COVID-19 is over, aircraft activity has been abnormal in May 2021 and from July to October 2023. We also find that the number of LTO and pollutant emissions of aircraft slightly differ at different hourly intervals, exhibiting a high rate in the daytime (06:00–23:00 UTC+8) and a low rate in the nighttime (00:00–05:00 UTC+8), with the minimum at 04:00 UTC+8 and the maximum at 13:00 UTC+8. In 2023, the aircraft activity was significantly better during the daytime (95.6 % of 2019 in LTO cycle) than in the nighttime (61.5 % of 2019 in LTO cycle).

In the LTO cycle, the horizontal distribution characteristics of aircraft emissions are dispersed along the runway, and the vertical distribution characteristics decrease as altitude increases. We find that aircraft emissions during the LTO cycle are so widely distributed around the airport that even a rough grid (e.g., 0.36° × 0.36°) cannot fully represent them. The elongated distribution characteristics of aircraft emissions indicate that evenly allocating emissions around the airport causes significant uncertainty. Due to variations in the MLH, the height at which aircraft emit pollutants during LTO can reach up to 4544 m above the maximum altitude of 915 m set by the ICAO. The NO_x emissions above the 915 m vary by month, accounting for an average of 24.6 % (9.9 %–37.5 %) in the LTO cycle.

Our 4D-LTO emission inventory dataset reflects the actual spatial and temporal and can be used to accurately assess the air quality impact of aircraft in the LTO cycle. This dataset and our methodology play a vital role in an in-depth study of temporal and spatial variations in aircraft emissions and their health and environmental impact. By conducting an in-depth analysis of our refined dataset, we can quantify the aviation industry's contribution to climate change and explore potential emission reduction pathways. Furthermore, by adjustments to accommodate regional differences, e.g., operational activity data, airport-specific emission factors, and airport-specific flight trajectory datasets, our methodology possesses broad applicability and flexibility. The application of our methodology to other regions is fundamental to formulating effective strategies and policies to achieve global aviation emission reduction targets.