# Development of Fast-Time Stochastic Airport Ground and Runway Simulation Model and Its Traffic Analysis

Authors: Ryota Mori
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Airport congestion, in particular congestion of departure aircraft, has already been discussed by other researches. Most solutions, though, fail to account for uncertainties. Since it is difficult to remove uncertainties of the operations in the real world, a strategy should be developed assuming such uncertainties exist. Therefore, this research develops a fast-time stochastic simulation model used to validate various methods in order to decrease airport congestion level under existing uncertainties. The surface movement data is analyzed first, and the uncertainty level is obtained. Next, based on the result of data analysis, the stochastic simulation model is developed. The model is validated statistically and the characteristics of airport operation under existing uncertainties are investigated.


---

## Body

## 1. Introduction

Airport ground congestions are becoming a critical problem at many airports in the world. Since the bottleneck of airport operations exists on the runway, there are long waiting queues of aircraft both on the ground taxiway and in the air, which increase the fuel burn and emissions. Arrival aircraft are often considered in research targeting airport congestions decrease because any additional flight time obviously requires extra fuel. Even if not so apparent, departure aircraft burn sufficient amount of fuel during taxiing, too, so departure queue management can help to reduce fuel burn. Although there are some researches regarding the taxi-out time saving of the departure aircraft [1–5], these researches focus on taxi-out time saving and do not investigate its negative effect. One such possible negative effect caused by departure queue management is take-off delay.The main reason for increased taxi-out time is that many departure aircraft wait in a queue before the runway due to runway congestion, so taxi-out time reduction is achieved by shifting the pushback time later intentionally. If the aircraft waits in the spot instead of waiting in a queue near the runway, the aircraft can turn its engines off and therefore save fuel. However, if this shift is too large, the aircraft cannot take off at the expected time. If all airport operations were estimated without errors, the reduction of taxi-out time would be maximized without imposing any take-off time delay, but this is impossible due to various uncertainties. Even if large margins are set to absorb uncertainties, the expected delay will be close to 0, but not definitely 0. Besides, setting of a large margin leads to decrease the reduction of taxi-out time as well. Therefore, to evaluate uncertainty effects, stochastic simulation model is necessary. The main focus of this paper is the development of such a stochastic simulation model.There are numerous airport simulation models proposed by many researchers already. However, this paper focuses on a stochastic model, which should also be appropriate to run a simulation fast enough. Most existing simulation models account for detailed aircraft movement but are also deterministic and slow, thus not suitable for the purpose of this paper [6–13]. Although some airport models consider uncertainty effect, such as the variance of taxiing speed or the take-off separation [14, 15], there are few researches considering uncertainties. In addition, uncertainty parameters are usually obtained via actual airport operation data, but only specific parameters are used in the simulation, and no stochastic simulation model has been verified whether it accurately models the airport operation. For example, even when the variation of taxiing speed is assumed, the model does not estimate the take-off time well unless the model handles variation of take-off separation. Therefore, this paper aims at developing a sufficiently accurate stochastic airport simulation model to account for the uncertainty effect. The stochastic parameters are obtained based on the actual operation data in each phase of the aircraft movement, and, using these parameters, the taxi-out of each departure aircraft is simulated. The simulated take-off time is then compared to the actual take-off time. Using the developed simulation model, the characteristics of the airport operation are also investigated, and the importance of the uncertainty is revealed.


## 2. Stochastic Airport Simulation Model at Tokyo International Airport

Tokyo International Airport is the target airport of this research. First, the airport operation is briefly explained in Section2.1, and a stochastic simulation model is developed. To conduct sufficient number of simulations to account for uncertainty effects, a single simulation run time should be short enough, so the model itself should be simplified. Error distribution models for stochastic components are introduced in Section 2.2, and taxi-out time of departure aircraft and taxi-in time of arrival aircraft are stochastically modeled in Section 2.3. Additional simulation constraints are explained in Section 2.4, and the take-off separation is stochastically modeled in Section 2.5. Finally, the simulation flow is shown and the model limitations are described in Section 2.6.


### 2.1. Tokyo International Airport and Its Runway Operation

Tokyo International Airport is the busiest airport in Japan and is mostly used for domestic flights. In 2010, the traffic volume was 303,000 flights per year and increases to 447,000 flights per year in 2014 with the opening of the new runway. Figure1 shows the airport map and typical operation under north wind. There are four runways at the airport. A runway is used for arrival and D runway is used for departure, but C runway is shared by both departure and arrival aircraft. Due to the runway location, aircraft departing from D runway cannot take off while a landing aircraft is approaching C runway. B runway is usually not used under north wind.Figure 1
The airport map and runway operation under north wind.In order to model taxiing correctly, knowledge of the procedures preceding take-off is necessary. The air traffic control (ATC) flow of departure aircraft is summarized in Figure2. First, about 5 minutes before the aircraft is ready for starting the engines, the pilot calls clearance delivery. If the flight plan is approved, the pilot will get a departure clearance from ATC. When the aircraft is ready for block-off, the pilot requests pushback to ATC. Once the pilot gets pushback approval, the aircraft starts pushback. During or after the pushback, the pilot requests taxiing to the runway. If the taxiing is approved and the aircraft is ready for taxiing, the aircraft will start taxiing. When the aircraft approaches the runway, the pilot requests runway clearance. Only after the runway clearance is approved, the aircraft can take off.Figure 2
Flow of departure.Here, several variables are defined. The time when the pilot starts pushback is AOBT (actual off-block time), and actual take-off time is defined as ATOT. The difference between ATOT and AOBT is defined as AXOT (actual taxi-out time). Airport operation is not completely deterministic, and there are probabilistic factors. Therefore, considering certain uncertainty, AXOT is determined probabilistically in the simulation. The distribution of AXOT is examined in Section2.3.As for arrival aircraft, the aircraft lands at ALDT (actual landing time). After landing, the aircraft go taxiing to the spot and block in at AIBT (actual in-block time). The duration of taxiing (AIBT-ALDT) is defined as AXIT (actual taxi-in time).This time, the data used to determine the simulation parameters in this research are obtained based on the smoothened airport surface movement data for 20 days between 2012 and 2014 (called Day 1 to Day 20), when north wind operation was conducted throughout a day.


### 2.2. Error Distribution Model

To consider the error factor, several error distribution models are applied. In this research, normal distribution and Erlang distribution are used. Normal distribution, also known as Gaussian distribution, is a symmetric distribution. Detailed explanation is not given here, but it has two parameters: averageμ and standard deviation (SD) σ. The probability density function of normal distribution is given by the following equation:(1)Nx;μ,σ=12πσexp-x-μ22σ2.Furthermore, to account for error’s asymmetry, Erlang distribution is also introduced. This distribution is asymmetric and is defined only whenx is greater than 0. This function is often used in the field of stochastic processes. There are two parameters: average μ and shape n (positive integer). The probability density function of Erlang distribution is given by the following equation. This distribution approaches the normal distribution as n increases:(2)Ex;μ,n=nnxn-1e-nx/μμnn-1!x>0.


### 2.3. Distribution of Taxi-Out Time (AXOT) and Taxi-In Time (AXIT)

Taxi-out time is defined as the time between pushback start and take-off. To determine the taxi-out time, it is divided into several stages, and the duration of each stage is determined. Figure3 shows the taxi-out flow. First, the aircraft has to complete pushback, defined as “pushback time” (Δtpushback). Next, the aircraft has to be released from the pushback truck and prepare for taxiing, defined as “preparation time” (Δtprepare). Then, the aircraft goes taxiing to the runway, defined as “taxiing time.” If there is a queue before the runway, the aircraft will need to wait extra. The minimum time which an aircraft needs to cover the distance between its spot and the runway even when no congestion is observed is defined as “minimum taxiing time” (Δttaxi). The difference between total taxiing time and minimum taxiing time is defined as “additional waiting time” (Δtwait). Finally, taxi-out time (AXOT) is calculated by the following equation. Each time is estimated based on the data available:(3)AXOT=Δtpushback+Δtprepare+Δttaxi+Δtwait.Here, some variables are defined. PTOT (earliest possible take-off time) is defined as the time when the aircraft is ready for take-off assuming there is no congestion. PTOT is calculated by the following equation:(4)PTOT=ATOT-Δtwait.Figure 3
Taxi-out flow.PXOT (earliest possible taxi-out time) is also defined in the same way by the following equation:(5)PXOT=AXOT-Δtwait=Δtpushback+Δtprepare+Δttaxi.As for the pushback time (Δtpushback), it is expected to depend on the pushback distance. The pushback distance and thus the pushback time are usually determined by the spot position. Another variable defined in this research is the preparation time (Δtprepare). It includes the time for the pilot to receive a taxiing clearance and the time needed to release the pushback truck from the aircraft. This preparation time is assumed to be the same in all situations, so the sum of the pushback time and the preparation time (Δtpushback+Δtprepare) is assumed to depend on the spot position, and this variable (called setup time) is used for the data analysis. Figure 4 shows the average and standard deviation of setup time in each spot position. As seen in the figure, the average of setup time varies with spot position, but the standard deviation of setup time does not change significantly with spot positions. Therefore, the setup time is assumed to depend on the spot position only. Here, the difference between the average setup time and the actual setup time is denoted by the residual of the setup time, and it is fitted by the combination of normal distribution and Erlang distribution as shown in Figure 5. The probability density function of setup time is calculated by the following function:(6)0.0959NΔtpushback+Δtprepare-t-setupspot;46.07,78.72+0.9041EΔtpushback+Δtprepare-t-setupspotktpushback+117.0;122.0,13,where t-setup(spot) is the average setup time in each spot obtained by data.Figure 4
Average and standard deviation of setup time in each spot position. (Spot number does not correspond to the actual spot number.)Figure 5
Residual of the setup time and fitted function.The next element needed to be defined in order to determine the earliest possible taxi-out time is the taxiing time (Δttaxi). The taxiing time is basically related to taxiing distance, so it is modeled with the parameter of the taxiing distance. Of course, the taxiing time varies with each pilot, and some aircraft go faster taxiing while others go slower. In addition, if two aircraft are conflicted along the taxiing route, additional taxiing time is required. Here, these effects are included in uncertainty. When obtaining the data of taxiing time, only the value of the sum of taxiing time and waiting time (Δttaxi+Δtwait) is available, because the data of start taxiing time and take-off time is obtained based on the surface movement data. Here, only the data where the aircraft goes taxiing smoothly and is not stuck in the waiting queue are used to calculate the taxiing time, because the waiting time is assumed to be zero for such aircraft.Figure6 shows the relationship between taxiing distance and taxiing time for noncongested aircraft departing from D runway. Note that the red points indicate the aircraft whose spot is at the international terminal. The figure shows that there is a correlation between taxiing distance and taxiing time, and the residual tends to increase with the taxiing distance. Therefore, the distribution of the normalized residual (residual per 1 km taxiing distance) is fitted by an Erlang distribution. In addition, the aircraft from the international terminal tends to have longer taxiing time, but that is because these aircraft have to cross A runway along the taxiing route and often wait due to passing landing aircraft. Therefore, it is assumed that these aircraft require additional taxiing time. Finally, the probability density function of taxiing time is calculated based on the following equations:(7)CRWY:EΔttaxi-0.1090dtaxi-68.03-Δx;37.3dtaxi1000,5,(8)DRWY:EΔttaxi-0.1002dtaxi-71.60-Δx;75.2dtaxi1000,30,(9)Δx=120,spotisintheinternationalterminal.0,otherwise.Figure 6
Taxiing time versus taxiing distance of aircraft departing from runway D.Taxiing route and route structure differ between C and D runways, so minimum taxiing time is calculated in a different manner. Once the spot position and the departure runway are determined, the taxiing distance and time are easily obtained. As for the arrival aircraft, the result of RWY C is well fitted, so the taxiing time is calculated based on (7).Finally, the additional waiting time (Δtwait) is considered. The additional waiting time is caused by the waiting queue at the runway, so it is strongly affected by take-off separation. Since only one aircraft can use the runway at the same time, a minimum separation (called take-off separation) is set. The take-off is usually operated based on first-come-first-served policy. If many aircraft come to the runway at the same time, a departure queue is made and the aircraft has to wait before the runway. In addition, the departure and arrival traffic are mutually dependent due to the arrangement of the runways at this airport, so the runway interaction should also be considered. The runway interaction is explained in Section 2.4. The take-off separation is affected by many parameters, so it is explained in Section 2.5.As for taxi-in time of arrival aircraft, the arrival aircraft only goes taxiing to the spot. If the spot is not occupied by other aircraft, the aircraft can block in. Therefore, uncertainty is found only in the taxiing phase. As mentioned before, the distribution of taxiing of arrival aircraft almost follows the one of departure aircraft on C runway, so (7) is used to estimate the taxiing time from the runway to the spot. The spot occupancy problem is described in Section 2.4.


### 2.4. Constraints at the Airport

In Section2.3, the normal operation of departure and arrival aircraft was explained. However, there are many constraints at the airport, such as (1) take-off separation, (2) mutual interaction between runways, (3) the spot occupancy problem, and (4) conflict of two aircraft on the taxiway. The run time of the simulation should be small, so this time only the constraints (1), (2), and (3) are considered, and the constraint (4) is not explicitly considered and is assumed to be included in “uncertainty.” The constraint (1) take-off separation will be explained in Section 2.5, so here the constraints (2) and (3) are explained.As for mutual interaction between runways, as shown in Figure1, the runway interaction is observed between take-off and landing aircraft on C runway and between take-off aircraft on D runway and landing aircraft on C runway. In addition, due to the departure and arrival route structure, C runway traffic and D runway traffic are mutually affected. However, this effect is complicated and its influence is relatively small so it is not considered here. Considering the runway operation, the landing time of arrival aircraft is currently not controlled, and the take-off time is controlled in accordance with the landing aircraft. Therefore, it is reasonable that the “no take-off time” relative to the landing time on C runway is set. According to the data analysis, the following constraints are set as shown in Table 1.Table 1
Prohibited departure time relative to the landing time on C runway.
Prohibited start
Prohibited end
Departure from C runway
−65 s
+85 s
Departure from D runway
−80 s
+0 sAs for the spot occupancy problem, obviously a single spot contains a single aircraft only. If the spot is already occupied by another aircraft, the arrival aircraft cannot block in until the aircraft leaves the spot, which is implemented in the simulation. In addition, paths near the spots are often shared among several spots, so arrival aircraft sometimes cannot get into the spot during the pushback of departure aircraft from a nearby spot. Therefore, based on the surface movement data, the conflict path of pushback and block-in is investigated in each spot, and the conflict effect between block-in path and pushback path is also implemented in the simulation.


### 2.5. Distribution of Take-Off Separation

The take-off separation is the key to determine the waiting time of departure aircraft. The take-off procedure is operated by a pilot, so the take-off separation should include uncertainty effect. However, the take-off separation is also affected by the following two factors: wake turbulence and weather condition.As for the wake turbulence, ICAO determines the minimum take-off separation based on the aircraft size of the current aircraft and the aircraft ahead [16]. There are four categories of aircraft: super, heavy, medium, and light (denoted by “S,” “H,” “M,” and “L”). Note that super and light aircraft are not operated at Tokyo International Airport, so only heavy and medium aircraft are considered in this research.Table2 shows the average take-off separation in each combination of wake turbulence category. Note that the standard deviation of take-off separation is about 20 s for each category. This table shows that the take-off separation differs by the wake turbulence category, and it is reflected in the calculation. The take-off separation between M-H and M-M is almost the same, so it is treated as the same value (90.75 s). From now, the take-off separation of “M-H” and “M-M” is treated as the nominal take-off separation, and the separation of “H-H” and “H-M” is reduced by the difference of the average to fit the nominal separation. The uncertainty of take-off separation is assumed to be the same for all wake turbulence categories.Table 2
Average take-off separation in each combination of wake turbulence category[s] (data size = 4144).
Preceding/following
Heavy
Medium
Heavy
99.92
106.06
Medium
90.28
91.31Regarding weather conditions, it is said that wind and visibility affect the take-off separation. According to the data analysis, only visibility statistically affects the take-off separation, and here only the visibility effect is explained. The visibility information is provided by METAR (METeological Airport Report) usually every thirty minutes. If the visibility is more than 10 km, it is recorded as 9999 m. Figure7 shows the relationship between the visibility and take-off separation (wake turbulence effect is already considered). Note that 80% of data is obtained when the visibility is more than 10 km, so low visibility data is relatively less. As shown in the figure, the take-off separation increases with smaller visibility, but it jumps up around 4000 m of visibility. Therefore, this is modeled by sigmoid function and linear regression as shown in the following expression:(10)95.64-0.0005356v+7.2351+exp0.0604v-4000,where v is the visibility in m. When the visibility of 10 km is treated as the nominal take-off separation, the separation data of another visibility is reduced to fit the nominal separation.Figure 7
Take-off separation versus visibility.Now, the nominal take-off separation is defined as the data where the visibility is 10 km and “M-H” or “M-M” of wake turbulence category is applied. Even if these effects are considered, there is still a large residual, which is modeled by the combination of normal distribution and Erlang distribution. Figure8 shows the distribution of obtained take-off separation and probability density function of the fitted distribution of the nominal take-off separation (tsep_nom) is described by the following equation:(11)0.890Etsep_nom-13.56;72.47,30+0.110Ntsep_nom;122.47,30.03.Figure 8
Distribution of take-off separation.The take-off separation should include the wake turbulence effect and visibility effect. Finally, the take-off separation (tsep) is calculated in the following equation:(12)tsep=tsep_nom-0.0005356v+7.2351+exp0.0604v-4000+5.35+Δtwake,Δtwake=9.17,H-H15.31,H-M0,M-HorM-M.


### 2.6. Flow of the Simulation and Limitations of the Simulation Model

In order to conduct a simulation, initial conditions must be specified. As for the departure aircraft, once AOBT is obtained, the PTOT can be calculated considering the uncertainty explained before. In the simulation, each departure aircraft is assumed to come to the runway at PTOT, and ATOT is decided based on the runway status and first-come-first-served basis. For arrival aircraft, once ALDT is obtained, AIBT is obtained based on (7) and taxiway status. The flow of the calculation of each variable is summarized in Figure 9.Figure 9
Simulation flow.In order to simulate the actual airport operation, AOBT and ALDT are set the same as the data obtained on each day. Furthermore, the spot position, the taxiing distance, the wake turbulence category, and the departure/arrival runway of each aircraft are also set based on each day’s data. These data are called each day’s scenario data. This scenario data includes each day’s traffic volume and the distribution of traffic. Even if the traffic volume is the same, the congestion level at the airport might differ between days. By using the scenario data, the daily fluctuation can also be investigated.In order to appropriately evaluate the results, the limitations of the simulation model should be carefully considered. First, this model does not take into account the conflict on taxiway between aircraft in any areas apart from the spot areas. It is assumed that the additional taxiing time by the conflict includes uncertainty. Second, in the real world, there are cases when the runway operation does not necessarily follow the first-come-first-served basis. The air traffic controllers are in charge of the operations, so the take-off sequence is sometimes changed. In addition, EDCT (expected departure clearance time) is sometimes set to the departure aircraft. This is the take-off time restriction due to the congestion in airspace or destination airport, and the aircraft cannot take off before EDCT, which is not considered in the simulation. EDCT also changes the departure sequence. Finally, nonstandard operation might be included in the data. It is confirmed in advance that the data do not include the long runway close, but short runway close might be included in data, which is difficult to exclude. Considering these limitations, the simulation accuracy will be evaluated.


## 2.1. Tokyo International Airport and Its Runway Operation

Tokyo International Airport is the busiest airport in Japan and is mostly used for domestic flights. In 2010, the traffic volume was 303,000 flights per year and increases to 447,000 flights per year in 2014 with the opening of the new runway. Figure1 shows the airport map and typical operation under north wind. There are four runways at the airport. A runway is used for arrival and D runway is used for departure, but C runway is shared by both departure and arrival aircraft. Due to the runway location, aircraft departing from D runway cannot take off while a landing aircraft is approaching C runway. B runway is usually not used under north wind.Figure 1
The airport map and runway operation under north wind.In order to model taxiing correctly, knowledge of the procedures preceding take-off is necessary. The air traffic control (ATC) flow of departure aircraft is summarized in Figure2. First, about 5 minutes before the aircraft is ready for starting the engines, the pilot calls clearance delivery. If the flight plan is approved, the pilot will get a departure clearance from ATC. When the aircraft is ready for block-off, the pilot requests pushback to ATC. Once the pilot gets pushback approval, the aircraft starts pushback. During or after the pushback, the pilot requests taxiing to the runway. If the taxiing is approved and the aircraft is ready for taxiing, the aircraft will start taxiing. When the aircraft approaches the runway, the pilot requests runway clearance. Only after the runway clearance is approved, the aircraft can take off.Figure 2
Flow of departure.Here, several variables are defined. The time when the pilot starts pushback is AOBT (actual off-block time), and actual take-off time is defined as ATOT. The difference between ATOT and AOBT is defined as AXOT (actual taxi-out time). Airport operation is not completely deterministic, and there are probabilistic factors. Therefore, considering certain uncertainty, AXOT is determined probabilistically in the simulation. The distribution of AXOT is examined in Section2.3.As for arrival aircraft, the aircraft lands at ALDT (actual landing time). After landing, the aircraft go taxiing to the spot and block in at AIBT (actual in-block time). The duration of taxiing (AIBT-ALDT) is defined as AXIT (actual taxi-in time).This time, the data used to determine the simulation parameters in this research are obtained based on the smoothened airport surface movement data for 20 days between 2012 and 2014 (called Day 1 to Day 20), when north wind operation was conducted throughout a day.


## 2.2. Error Distribution Model

To consider the error factor, several error distribution models are applied. In this research, normal distribution and Erlang distribution are used. Normal distribution, also known as Gaussian distribution, is a symmetric distribution. Detailed explanation is not given here, but it has two parameters: averageμ and standard deviation (SD) σ. The probability density function of normal distribution is given by the following equation:(1)Nx;μ,σ=12πσexp-x-μ22σ2.Furthermore, to account for error’s asymmetry, Erlang distribution is also introduced. This distribution is asymmetric and is defined only whenx is greater than 0. This function is often used in the field of stochastic processes. There are two parameters: average μ and shape n (positive integer). The probability density function of Erlang distribution is given by the following equation. This distribution approaches the normal distribution as n increases:(2)Ex;μ,n=nnxn-1e-nx/μμnn-1!x>0.


## 2.3. Distribution of Taxi-Out Time (AXOT) and Taxi-In Time (AXIT)

Taxi-out time is defined as the time between pushback start and take-off. To determine the taxi-out time, it is divided into several stages, and the duration of each stage is determined. Figure3 shows the taxi-out flow. First, the aircraft has to complete pushback, defined as “pushback time” (Δtpushback). Next, the aircraft has to be released from the pushback truck and prepare for taxiing, defined as “preparation time” (Δtprepare). Then, the aircraft goes taxiing to the runway, defined as “taxiing time.” If there is a queue before the runway, the aircraft will need to wait extra. The minimum time which an aircraft needs to cover the distance between its spot and the runway even when no congestion is observed is defined as “minimum taxiing time” (Δttaxi). The difference between total taxiing time and minimum taxiing time is defined as “additional waiting time” (Δtwait). Finally, taxi-out time (AXOT) is calculated by the following equation. Each time is estimated based on the data available:(3)AXOT=Δtpushback+Δtprepare+Δttaxi+Δtwait.Here, some variables are defined. PTOT (earliest possible take-off time) is defined as the time when the aircraft is ready for take-off assuming there is no congestion. PTOT is calculated by the following equation:(4)PTOT=ATOT-Δtwait.Figure 3
Taxi-out flow.PXOT (earliest possible taxi-out time) is also defined in the same way by the following equation:(5)PXOT=AXOT-Δtwait=Δtpushback+Δtprepare+Δttaxi.As for the pushback time (Δtpushback), it is expected to depend on the pushback distance. The pushback distance and thus the pushback time are usually determined by the spot position. Another variable defined in this research is the preparation time (Δtprepare). It includes the time for the pilot to receive a taxiing clearance and the time needed to release the pushback truck from the aircraft. This preparation time is assumed to be the same in all situations, so the sum of the pushback time and the preparation time (Δtpushback+Δtprepare) is assumed to depend on the spot position, and this variable (called setup time) is used for the data analysis. Figure 4 shows the average and standard deviation of setup time in each spot position. As seen in the figure, the average of setup time varies with spot position, but the standard deviation of setup time does not change significantly with spot positions. Therefore, the setup time is assumed to depend on the spot position only. Here, the difference between the average setup time and the actual setup time is denoted by the residual of the setup time, and it is fitted by the combination of normal distribution and Erlang distribution as shown in Figure 5. The probability density function of setup time is calculated by the following function:(6)0.0959NΔtpushback+Δtprepare-t-setupspot;46.07,78.72+0.9041EΔtpushback+Δtprepare-t-setupspotktpushback+117.0;122.0,13,where t-setup(spot) is the average setup time in each spot obtained by data.Figure 4
Average and standard deviation of setup time in each spot position. (Spot number does not correspond to the actual spot number.)Figure 5
Residual of the setup time and fitted function.The next element needed to be defined in order to determine the earliest possible taxi-out time is the taxiing time (Δttaxi). The taxiing time is basically related to taxiing distance, so it is modeled with the parameter of the taxiing distance. Of course, the taxiing time varies with each pilot, and some aircraft go faster taxiing while others go slower. In addition, if two aircraft are conflicted along the taxiing route, additional taxiing time is required. Here, these effects are included in uncertainty. When obtaining the data of taxiing time, only the value of the sum of taxiing time and waiting time (Δttaxi+Δtwait) is available, because the data of start taxiing time and take-off time is obtained based on the surface movement data. Here, only the data where the aircraft goes taxiing smoothly and is not stuck in the waiting queue are used to calculate the taxiing time, because the waiting time is assumed to be zero for such aircraft.Figure6 shows the relationship between taxiing distance and taxiing time for noncongested aircraft departing from D runway. Note that the red points indicate the aircraft whose spot is at the international terminal. The figure shows that there is a correlation between taxiing distance and taxiing time, and the residual tends to increase with the taxiing distance. Therefore, the distribution of the normalized residual (residual per 1 km taxiing distance) is fitted by an Erlang distribution. In addition, the aircraft from the international terminal tends to have longer taxiing time, but that is because these aircraft have to cross A runway along the taxiing route and often wait due to passing landing aircraft. Therefore, it is assumed that these aircraft require additional taxiing time. Finally, the probability density function of taxiing time is calculated based on the following equations:(7)CRWY:EΔttaxi-0.1090dtaxi-68.03-Δx;37.3dtaxi1000,5,(8)DRWY:EΔttaxi-0.1002dtaxi-71.60-Δx;75.2dtaxi1000,30,(9)Δx=120,spotisintheinternationalterminal.0,otherwise.Figure 6
Taxiing time versus taxiing distance of aircraft departing from runway D.Taxiing route and route structure differ between C and D runways, so minimum taxiing time is calculated in a different manner. Once the spot position and the departure runway are determined, the taxiing distance and time are easily obtained. As for the arrival aircraft, the result of RWY C is well fitted, so the taxiing time is calculated based on (7).Finally, the additional waiting time (Δtwait) is considered. The additional waiting time is caused by the waiting queue at the runway, so it is strongly affected by take-off separation. Since only one aircraft can use the runway at the same time, a minimum separation (called take-off separation) is set. The take-off is usually operated based on first-come-first-served policy. If many aircraft come to the runway at the same time, a departure queue is made and the aircraft has to wait before the runway. In addition, the departure and arrival traffic are mutually dependent due to the arrangement of the runways at this airport, so the runway interaction should also be considered. The runway interaction is explained in Section 2.4. The take-off separation is affected by many parameters, so it is explained in Section 2.5.As for taxi-in time of arrival aircraft, the arrival aircraft only goes taxiing to the spot. If the spot is not occupied by other aircraft, the aircraft can block in. Therefore, uncertainty is found only in the taxiing phase. As mentioned before, the distribution of taxiing of arrival aircraft almost follows the one of departure aircraft on C runway, so (7) is used to estimate the taxiing time from the runway to the spot. The spot occupancy problem is described in Section 2.4.


## 2.4. Constraints at the Airport

In Section2.3, the normal operation of departure and arrival aircraft was explained. However, there are many constraints at the airport, such as (1) take-off separation, (2) mutual interaction between runways, (3) the spot occupancy problem, and (4) conflict of two aircraft on the taxiway. The run time of the simulation should be small, so this time only the constraints (1), (2), and (3) are considered, and the constraint (4) is not explicitly considered and is assumed to be included in “uncertainty.” The constraint (1) take-off separation will be explained in Section 2.5, so here the constraints (2) and (3) are explained.As for mutual interaction between runways, as shown in Figure1, the runway interaction is observed between take-off and landing aircraft on C runway and between take-off aircraft on D runway and landing aircraft on C runway. In addition, due to the departure and arrival route structure, C runway traffic and D runway traffic are mutually affected. However, this effect is complicated and its influence is relatively small so it is not considered here. Considering the runway operation, the landing time of arrival aircraft is currently not controlled, and the take-off time is controlled in accordance with the landing aircraft. Therefore, it is reasonable that the “no take-off time” relative to the landing time on C runway is set. According to the data analysis, the following constraints are set as shown in Table 1.Table 1
Prohibited departure time relative to the landing time on C runway.
Prohibited start
Prohibited end
Departure from C runway
−65 s
+85 s
Departure from D runway
−80 s
+0 sAs for the spot occupancy problem, obviously a single spot contains a single aircraft only. If the spot is already occupied by another aircraft, the arrival aircraft cannot block in until the aircraft leaves the spot, which is implemented in the simulation. In addition, paths near the spots are often shared among several spots, so arrival aircraft sometimes cannot get into the spot during the pushback of departure aircraft from a nearby spot. Therefore, based on the surface movement data, the conflict path of pushback and block-in is investigated in each spot, and the conflict effect between block-in path and pushback path is also implemented in the simulation.


## 2.5. Distribution of Take-Off Separation

The take-off separation is the key to determine the waiting time of departure aircraft. The take-off procedure is operated by a pilot, so the take-off separation should include uncertainty effect. However, the take-off separation is also affected by the following two factors: wake turbulence and weather condition.As for the wake turbulence, ICAO determines the minimum take-off separation based on the aircraft size of the current aircraft and the aircraft ahead [16]. There are four categories of aircraft: super, heavy, medium, and light (denoted by “S,” “H,” “M,” and “L”). Note that super and light aircraft are not operated at Tokyo International Airport, so only heavy and medium aircraft are considered in this research.Table2 shows the average take-off separation in each combination of wake turbulence category. Note that the standard deviation of take-off separation is about 20 s for each category. This table shows that the take-off separation differs by the wake turbulence category, and it is reflected in the calculation. The take-off separation between M-H and M-M is almost the same, so it is treated as the same value (90.75 s). From now, the take-off separation of “M-H” and “M-M” is treated as the nominal take-off separation, and the separation of “H-H” and “H-M” is reduced by the difference of the average to fit the nominal separation. The uncertainty of take-off separation is assumed to be the same for all wake turbulence categories.Table 2
Average take-off separation in each combination of wake turbulence category[s] (data size = 4144).
Preceding/following
Heavy
Medium
Heavy
99.92
106.06
Medium
90.28
91.31Regarding weather conditions, it is said that wind and visibility affect the take-off separation. According to the data analysis, only visibility statistically affects the take-off separation, and here only the visibility effect is explained. The visibility information is provided by METAR (METeological Airport Report) usually every thirty minutes. If the visibility is more than 10 km, it is recorded as 9999 m. Figure7 shows the relationship between the visibility and take-off separation (wake turbulence effect is already considered). Note that 80% of data is obtained when the visibility is more than 10 km, so low visibility data is relatively less. As shown in the figure, the take-off separation increases with smaller visibility, but it jumps up around 4000 m of visibility. Therefore, this is modeled by sigmoid function and linear regression as shown in the following expression:(10)95.64-0.0005356v+7.2351+exp0.0604v-4000,where v is the visibility in m. When the visibility of 10 km is treated as the nominal take-off separation, the separation data of another visibility is reduced to fit the nominal separation.Figure 7
Take-off separation versus visibility.Now, the nominal take-off separation is defined as the data where the visibility is 10 km and “M-H” or “M-M” of wake turbulence category is applied. Even if these effects are considered, there is still a large residual, which is modeled by the combination of normal distribution and Erlang distribution. Figure8 shows the distribution of obtained take-off separation and probability density function of the fitted distribution of the nominal take-off separation (tsep_nom) is described by the following equation:(11)0.890Etsep_nom-13.56;72.47,30+0.110Ntsep_nom;122.47,30.03.Figure 8
Distribution of take-off separation.The take-off separation should include the wake turbulence effect and visibility effect. Finally, the take-off separation (tsep) is calculated in the following equation:(12)tsep=tsep_nom-0.0005356v+7.2351+exp0.0604v-4000+5.35+Δtwake,Δtwake=9.17,H-H15.31,H-M0,M-HorM-M.


## 2.6. Flow of the Simulation and Limitations of the Simulation Model

In order to conduct a simulation, initial conditions must be specified. As for the departure aircraft, once AOBT is obtained, the PTOT can be calculated considering the uncertainty explained before. In the simulation, each departure aircraft is assumed to come to the runway at PTOT, and ATOT is decided based on the runway status and first-come-first-served basis. For arrival aircraft, once ALDT is obtained, AIBT is obtained based on (7) and taxiway status. The flow of the calculation of each variable is summarized in Figure 9.Figure 9
Simulation flow.In order to simulate the actual airport operation, AOBT and ALDT are set the same as the data obtained on each day. Furthermore, the spot position, the taxiing distance, the wake turbulence category, and the departure/arrival runway of each aircraft are also set based on each day’s data. These data are called each day’s scenario data. This scenario data includes each day’s traffic volume and the distribution of traffic. Even if the traffic volume is the same, the congestion level at the airport might differ between days. By using the scenario data, the daily fluctuation can also be investigated.In order to appropriately evaluate the results, the limitations of the simulation model should be carefully considered. First, this model does not take into account the conflict on taxiway between aircraft in any areas apart from the spot areas. It is assumed that the additional taxiing time by the conflict includes uncertainty. Second, in the real world, there are cases when the runway operation does not necessarily follow the first-come-first-served basis. The air traffic controllers are in charge of the operations, so the take-off sequence is sometimes changed. In addition, EDCT (expected departure clearance time) is sometimes set to the departure aircraft. This is the take-off time restriction due to the congestion in airspace or destination airport, and the aircraft cannot take off before EDCT, which is not considered in the simulation. EDCT also changes the departure sequence. Finally, nonstandard operation might be included in the data. It is confirmed in advance that the data do not include the long runway close, but short runway close might be included in data, which is difficult to exclude. Considering these limitations, the simulation accuracy will be evaluated.


## 3. Verification of Proposed Simulation Model and Daily Data Analysis

### 3.1. Characteristics of Daily Data

Before evaluating the simulation model, the characteristics of the daily traffic are investigated first. The proposed simulation model is to be used to evaluate uncertainty effect, especially important in congestions, so the model performance in such cases is of the utmost importance. Therefore, first, the congestion level throughout a day at the airport is investigated.When the airport is not congested, the take-off is operated smoothly; that is, the waiting time of take-off aircraft is zero. Therefore, the congestion level at the airport closely relates to the waiting time of take-off aircraft. To investigate the actual congestion level at the airport, the actual waiting time of take-off aircraft is calculated. Since the actual waiting time is difficult to obtain directly, it is estimated based on the difference between AXOT and PXOT. AXOT can be obtained directly from the data. PXOT can be estimated as the nominal PXOT. PXOT is usually stochastically calculated as explained in the last section, but the nominal PXOT can be estimated if the uncertainty is assumed to be zero. Even if the uncertainty is zero, the nominal PXOT includes the taxiing distance effect and the spot position effect. In this way, the waiting time is estimated for each aircraft, and this waiting time based on actual data is called estimated waiting time. Figure10 shows the estimated waiting time of each aircraft throughout a day on Day 1. The waiting time is not distributed evenly throughout a day, because the traffic volume in each time range differs. Some data include negative waiting time, which occurs because the estimated waiting time is calculated based on the estimated PXOT. Now, the time is split into three hours each, and the congestion is considered in each time range. The traffic volume and scheduled traffic distribution are almost the same within the time range between days, so the waiting time can be compared in each time range.Figure 10
Estimated waiting time on Day 1 (ordered by actual take-off sequence).Figure11 shows the total estimated waiting time in each time range on each day. Table 3 shows the average of estimated total waiting time for 20 days and the average number of departure and arrival aircraft. According to the figure and the table, the largest total waiting time is observed at PM6–PM9, and the smallest waiting time is observed at PM3–PM6, but the daily fluctuation is also large. On the other hand, both departure and arrival traffic are the largest at AM9–AM12. Although more traffic potentially causes longer waiting time, there seems to be other factors as well. The waiting time increases if the runway waiting queue gets longer. The length of the runway queue is heavily affected by the traffic concentration. If many aircraft approach the runway at the same time, the runway queue becomes long and therefore the waiting time gets long and vice versa. Actually, at PM6–PM9, the departure traffic is not evenly distributed for these 3 hours, and it is concentrated between PM7 and PM8 and the waiting time is also long. Therefore, the effect of the traffic concentration seems to be more important than the traffic volume itself. Each day’s scenario fixes AOBT, which includes the actual distribution of off-block time, so the traffic concentration effect can also be evaluated in the simulation.Table 3
Average estimated waiting time and traffic volume in each time range for 20 days.
Time range
Estimated waiting time [minutes]
Number of departure aircraft
Number of arrival aircraft
AM6–AM9
174.5
105.5
44.0
AM9–AM12
178.3
105.8
104.4
AM12–PM3
158.3
99.6
98.8
PM3–PM6
133.4
93.8
102.1
PM6–PM9
281.9
98.4
104.2Figure 11
Total estimated waiting time in each time range on each day.In addition, the take-off separation can also affect the total waiting time. According to Figure7 and (11), the visibility affects the take-off separation. When the take-off separation is large, large total waiting time is expected. In Figure 11, the cases where the visibility is less than 5000 m are shown separately. When the visibility is less than 5000 m, large total waiting time is observed at all times. However, the cases where the large total waiting time is observed are not necessarily on a low visibility day. The visibility, that is, the take-off separation, is an important factor, but it is not the only one.


### 3.2. Validation of the Proposed Simulation Model

Next, the proposed simulation model is evaluated. The model considers the uncertainty effect, so it should be evaluated by a sufficient number of simulation runs. This time, 10,000 runs of simulations are conducted and the result is discussed. The program is made by C++ language and run with Intel Core i7-3770. It takes about 20 s to complete 10,000 runs of simulation on each day, and it is sufficiently fast to evaluate the uncertainty effect.When evaluating the simulation model, it is important to decide what is expected in the simulation model. The main purpose of the simulation model is to evaluate the uncertainty effect, and the modeling of the congestion phenomena is the most important. On the other hand, the runway operation is assumed to follow the first-come-first-served basis, so the take-off time of individual aircraft (i.e., waiting time of individual aircraft) might not be necessarily well modeled. Here, the total waiting time and the take-off time of each aircraft are used as an index of the simulation accuracy.First, the total waiting time is discussed. Figure12 shows the distribution of the total waiting time on Day 1 PM6–PM9 based on 10,000 times of simulation. The estimated waiting time based on the actual data is 316.5 minutes, and the peak of the distribution by the simulation and the estimated waiting time almost match. The estimated waiting time corresponds to 47.49 percentile of the distribution. However, only this result by itself does not indicate that the proposed simulation model works well. The actual operation can be considered a single sample set of simulations. If so, it is expected that the actual operation is the same as one of the simulations. To evaluate it, the percentile of the waiting time is used. If the uncertainty of the simulation is underestimated, the percentile will be observed often around 0 and 100. If the uncertainty is overestimated, the percentile will be observed often around 50. Therefore, if the uncertainty of the simulation is well modeled, the percentile of the estimated waiting time should be uniformly distributed between 0 and 100. As shown in Figure 11, a whole day is split into 5 time ranges for all 20 days; the percentile of the total estimated waiting time of 100 data is obtained. To investigate that the obtained percentile is evenly distributed, Q-Q plot is often used. Q-Q plot is a probability plot, where one axis shows the data quantile and the other axis shows the theoretical quantile. If the actual data completely follow the theoretical distribution, the plot is observed on the line “y=x.” Figure 13 shows the Q-Q plot of the total waiting time.Figure 12
The distribution of total waiting time based on 10,000 runs of simulation and the estimated waiting time on Day 1 PM6–PM9.Figure 13
Q-Q plot of the total estimated waiting time on each day and each time range.This figure shows that the theoretical quantile and data quantile match very well, which infers that the uncertainty is well modeled, and therefore the proposed simulation model works well. The uniformity of the distribution is also examined via a statistical test. This time, Kolmogorov-Smirnov test (KS test) [17] and Anderson-Darling test (AD test) [18] are used. Both tests are statistical tests to verify whether a given sample of data is made from a given probability distribution (this time, uniform distribution). Both tests show a P value, which is the probability of a test result being at least as extreme as the one that is actually obtained. A small P value means that the actual data is more extreme, and usually if the P value is less than 0.05, it is concluded that the actual data do not come from a given probability distribution. The main difference between KS test and AD test is that AD test weights much more on the tail probability. However, these tests are usually done to reject the null hypothesis, that is, to prove that the obtained data do not come from the given probability distribution. Therefore, even if the P value is greater than 0.05, it does not directly mean that the obtained data follow the given probability distribution, but such a result leaves open the probability that the obtained data do not come from the given probability distribution.The statistical test results are shown in Table4. The statistical test is also done with data of each time range only. According to the result, the P value of all data is much greater than 0.05 for both tests, and no P value being less than 0.05 is observed for any time range.Table 4
Statistical tests results.
P value
KS test
AD test
AM6–AM9
0.8247
0.8443
AM9–PM12
0.2684
0.1409
PM12–PM3
0.7513
0.6186
PM3–PM6
0.4026
0.0693
PM6–PM9
0.3282
0.5481
All
0.8211
0.4698Next, the take-off time is considered. The take-off time in the simulation is also distributed, and the percentile of the actual take-off time is obtained. However, as mentioned before, the take-off sequence is not modeled in the simulation, so the percentile of the take-off time might not be uniformly distributed. Figure14 shows the Q-Q plot of the take-off time of all aircraft for 20 days. In total, 10,226 departure aircraft data are used.Figure 14
Q-Q plot of take-off time.The data quantile slightly differs from the theoretical quantile around the edge of the quantile. KS test is also conducted, and theP value here is 1.385E-9, which means that the actual percentile is not uniformly distributed. However, as mentioned in the last section, this simulation model does not include all effects of airport operation, and it is expected that the percentile of take-off time is not completely uniformly distributed. However, according to Figure 14, the data quantile and theoretical quartile almost match each other, and also considering the result of total waiting time, it is concluded that the proposed simulation model is sufficiently accurate as an airport-runway statistical simulation model.


### 3.3. Uncertainty Effect Based on the Simulation

Since the simulation model is developed successfully, several aspects of uncertainty effects in airport operations are examined. The congestion in a runway queue is observed every day, but the congestion level changes every day. There are many possible factors related to the total waiting time, such as the number of departure aircraft, the initial condition (how much traffic is concentrated on a specific time slot, here also called scenario uncertainty), uncertainties of taxi-in, taxi-out, and take-off separation (called traffic uncertainty), and meteorological conditions. The effects of the number of departure aircraft and meteorological conditions can be associated with each day. The traffic uncertainty effect can be evaluated with a large number of simulations under the same initial condition. The scenario uncertainty can be investigated between days. It is interesting to understand how much each factor affects the congestion level.First, the variation of the congestion within a day is examined. An example is shown in Figure12 on Day PM6–PM9, which shows the distribution of total waiting time based on 10,000 simulations. According to additional calculations, the similar shape of distribution is obtained on other days and other time ranges. This variation of the total waiting time stems from traffic uncertainty effects. On this day, the estimated total waiting time based on the actual data is 316.5 minutes, but the distribution of the total waiting time between about 150 minutes and 600 minutes is observed. Even between 2.5 percentile and 97.5 percentile, the total waiting time varies between 211 and 472 minutes. This infers that the congestion level changes more than twice even if the scenario involves no uncertainty. In addition, the distribution of the total waiting time is not symmetric, and the large waiting time is observed slightly more often.Next, Figure15 shows the daily fluctuation of the total waiting time at PM6–PM9. The bar indicates the estimated total waiting time based on actual data (the same data shown in Figure 11), and the error bar indicates the 95% range of the total waiting time in the simulation. The red line shows the number of departure aircraft within this time range. On Day 15, visibility was low. From this figure, there are many interesting points found.Figure 15
Total estimated waiting time and 95% range of total simulated waiting time at PM6–PM9 on each day.First, the number of departure aircraft is almost constant in this time range, and only a small difference among the days is observed. However, like the result obtained in Table3, there is no clear relationship between the number of departure aircraft and the total waiting time. Second, even though the flight schedule is almost the same on each day in each time range, the range of the total waiting time in the simulation varies. The simulation scenario includes the initial condition of the off-block time, so it does not include the uncertainty between actual off-block time and scheduled off-block time, that is, scenario uncertainty. On Days 15 and 17, large estimated waiting time is observed, but the simulated total waiting time is also large. The low visibility seems to be a reason for large total waiting time on Day 15, but that is not the case on Day 17, maybe due to the traffic concentration. Therefore, on these days, the large total waiting time is caused not because of the traffic uncertainty effect, but because of either visibility or the scenario uncertainty, that is, the uncertainty of the off-block time compared to the scheduled time. If so, the total waiting time in the simulation varies very much with only traffic uncertainty, which means that the further large distribution will be obtained if the scenario uncertainty is considered. Therefore, the uncertainty affects the airport congestion very much, and it is almost impossible to estimate the airport congestion level on a specific day in advance, though the average congestion level can be estimated.However, this result does not directly mean that congestion cannot be relieved. If the scheduled departure or arrival time is optimally decided, the average congestion level might be reduced with keeping the traffic volume. The better spot allocation might also decrease the congestion level. However, since we understand that the uncertainty in airport traffic is significant, the uncertainty effect cannot be ignored when considering the airport operation. Otherwise, even if the taxi-out time is reduced, the uncertainty might cause the delay of take-off time or reduce the runway capacity.


## 3.1. Characteristics of Daily Data

Before evaluating the simulation model, the characteristics of the daily traffic are investigated first. The proposed simulation model is to be used to evaluate uncertainty effect, especially important in congestions, so the model performance in such cases is of the utmost importance. Therefore, first, the congestion level throughout a day at the airport is investigated.When the airport is not congested, the take-off is operated smoothly; that is, the waiting time of take-off aircraft is zero. Therefore, the congestion level at the airport closely relates to the waiting time of take-off aircraft. To investigate the actual congestion level at the airport, the actual waiting time of take-off aircraft is calculated. Since the actual waiting time is difficult to obtain directly, it is estimated based on the difference between AXOT and PXOT. AXOT can be obtained directly from the data. PXOT can be estimated as the nominal PXOT. PXOT is usually stochastically calculated as explained in the last section, but the nominal PXOT can be estimated if the uncertainty is assumed to be zero. Even if the uncertainty is zero, the nominal PXOT includes the taxiing distance effect and the spot position effect. In this way, the waiting time is estimated for each aircraft, and this waiting time based on actual data is called estimated waiting time. Figure10 shows the estimated waiting time of each aircraft throughout a day on Day 1. The waiting time is not distributed evenly throughout a day, because the traffic volume in each time range differs. Some data include negative waiting time, which occurs because the estimated waiting time is calculated based on the estimated PXOT. Now, the time is split into three hours each, and the congestion is considered in each time range. The traffic volume and scheduled traffic distribution are almost the same within the time range between days, so the waiting time can be compared in each time range.Figure 10
Estimated waiting time on Day 1 (ordered by actual take-off sequence).Figure11 shows the total estimated waiting time in each time range on each day. Table 3 shows the average of estimated total waiting time for 20 days and the average number of departure and arrival aircraft. According to the figure and the table, the largest total waiting time is observed at PM6–PM9, and the smallest waiting time is observed at PM3–PM6, but the daily fluctuation is also large. On the other hand, both departure and arrival traffic are the largest at AM9–AM12. Although more traffic potentially causes longer waiting time, there seems to be other factors as well. The waiting time increases if the runway waiting queue gets longer. The length of the runway queue is heavily affected by the traffic concentration. If many aircraft approach the runway at the same time, the runway queue becomes long and therefore the waiting time gets long and vice versa. Actually, at PM6–PM9, the departure traffic is not evenly distributed for these 3 hours, and it is concentrated between PM7 and PM8 and the waiting time is also long. Therefore, the effect of the traffic concentration seems to be more important than the traffic volume itself. Each day’s scenario fixes AOBT, which includes the actual distribution of off-block time, so the traffic concentration effect can also be evaluated in the simulation.Table 3
Average estimated waiting time and traffic volume in each time range for 20 days.
Time range
Estimated waiting time [minutes]
Number of departure aircraft
Number of arrival aircraft
AM6–AM9
174.5
105.5
44.0
AM9–AM12
178.3
105.8
104.4
AM12–PM3
158.3
99.6
98.8
PM3–PM6
133.4
93.8
102.1
PM6–PM9
281.9
98.4
104.2Figure 11
Total estimated waiting time in each time range on each day.In addition, the take-off separation can also affect the total waiting time. According to Figure7 and (11), the visibility affects the take-off separation. When the take-off separation is large, large total waiting time is expected. In Figure 11, the cases where the visibility is less than 5000 m are shown separately. When the visibility is less than 5000 m, large total waiting time is observed at all times. However, the cases where the large total waiting time is observed are not necessarily on a low visibility day. The visibility, that is, the take-off separation, is an important factor, but it is not the only one.


## 3.2. Validation of the Proposed Simulation Model

Next, the proposed simulation model is evaluated. The model considers the uncertainty effect, so it should be evaluated by a sufficient number of simulation runs. This time, 10,000 runs of simulations are conducted and the result is discussed. The program is made by C++ language and run with Intel Core i7-3770. It takes about 20 s to complete 10,000 runs of simulation on each day, and it is sufficiently fast to evaluate the uncertainty effect.When evaluating the simulation model, it is important to decide what is expected in the simulation model. The main purpose of the simulation model is to evaluate the uncertainty effect, and the modeling of the congestion phenomena is the most important. On the other hand, the runway operation is assumed to follow the first-come-first-served basis, so the take-off time of individual aircraft (i.e., waiting time of individual aircraft) might not be necessarily well modeled. Here, the total waiting time and the take-off time of each aircraft are used as an index of the simulation accuracy.First, the total waiting time is discussed. Figure12 shows the distribution of the total waiting time on Day 1 PM6–PM9 based on 10,000 times of simulation. The estimated waiting time based on the actual data is 316.5 minutes, and the peak of the distribution by the simulation and the estimated waiting time almost match. The estimated waiting time corresponds to 47.49 percentile of the distribution. However, only this result by itself does not indicate that the proposed simulation model works well. The actual operation can be considered a single sample set of simulations. If so, it is expected that the actual operation is the same as one of the simulations. To evaluate it, the percentile of the waiting time is used. If the uncertainty of the simulation is underestimated, the percentile will be observed often around 0 and 100. If the uncertainty is overestimated, the percentile will be observed often around 50. Therefore, if the uncertainty of the simulation is well modeled, the percentile of the estimated waiting time should be uniformly distributed between 0 and 100. As shown in Figure 11, a whole day is split into 5 time ranges for all 20 days; the percentile of the total estimated waiting time of 100 data is obtained. To investigate that the obtained percentile is evenly distributed, Q-Q plot is often used. Q-Q plot is a probability plot, where one axis shows the data quantile and the other axis shows the theoretical quantile. If the actual data completely follow the theoretical distribution, the plot is observed on the line “y=x.” Figure 13 shows the Q-Q plot of the total waiting time.Figure 12
The distribution of total waiting time based on 10,000 runs of simulation and the estimated waiting time on Day 1 PM6–PM9.Figure 13
Q-Q plot of the total estimated waiting time on each day and each time range.This figure shows that the theoretical quantile and data quantile match very well, which infers that the uncertainty is well modeled, and therefore the proposed simulation model works well. The uniformity of the distribution is also examined via a statistical test. This time, Kolmogorov-Smirnov test (KS test) [17] and Anderson-Darling test (AD test) [18] are used. Both tests are statistical tests to verify whether a given sample of data is made from a given probability distribution (this time, uniform distribution). Both tests show a P value, which is the probability of a test result being at least as extreme as the one that is actually obtained. A small P value means that the actual data is more extreme, and usually if the P value is less than 0.05, it is concluded that the actual data do not come from a given probability distribution. The main difference between KS test and AD test is that AD test weights much more on the tail probability. However, these tests are usually done to reject the null hypothesis, that is, to prove that the obtained data do not come from the given probability distribution. Therefore, even if the P value is greater than 0.05, it does not directly mean that the obtained data follow the given probability distribution, but such a result leaves open the probability that the obtained data do not come from the given probability distribution.The statistical test results are shown in Table4. The statistical test is also done with data of each time range only. According to the result, the P value of all data is much greater than 0.05 for both tests, and no P value being less than 0.05 is observed for any time range.Table 4
Statistical tests results.
P value
KS test
AD test
AM6–AM9
0.8247
0.8443
AM9–PM12
0.2684
0.1409
PM12–PM3
0.7513
0.6186
PM3–PM6
0.4026
0.0693
PM6–PM9
0.3282
0.5481
All
0.8211
0.4698Next, the take-off time is considered. The take-off time in the simulation is also distributed, and the percentile of the actual take-off time is obtained. However, as mentioned before, the take-off sequence is not modeled in the simulation, so the percentile of the take-off time might not be uniformly distributed. Figure14 shows the Q-Q plot of the take-off time of all aircraft for 20 days. In total, 10,226 departure aircraft data are used.Figure 14
Q-Q plot of take-off time.The data quantile slightly differs from the theoretical quantile around the edge of the quantile. KS test is also conducted, and theP value here is 1.385E-9, which means that the actual percentile is not uniformly distributed. However, as mentioned in the last section, this simulation model does not include all effects of airport operation, and it is expected that the percentile of take-off time is not completely uniformly distributed. However, according to Figure 14, the data quantile and theoretical quartile almost match each other, and also considering the result of total waiting time, it is concluded that the proposed simulation model is sufficiently accurate as an airport-runway statistical simulation model.


## 3.3. Uncertainty Effect Based on the Simulation

Since the simulation model is developed successfully, several aspects of uncertainty effects in airport operations are examined. The congestion in a runway queue is observed every day, but the congestion level changes every day. There are many possible factors related to the total waiting time, such as the number of departure aircraft, the initial condition (how much traffic is concentrated on a specific time slot, here also called scenario uncertainty), uncertainties of taxi-in, taxi-out, and take-off separation (called traffic uncertainty), and meteorological conditions. The effects of the number of departure aircraft and meteorological conditions can be associated with each day. The traffic uncertainty effect can be evaluated with a large number of simulations under the same initial condition. The scenario uncertainty can be investigated between days. It is interesting to understand how much each factor affects the congestion level.First, the variation of the congestion within a day is examined. An example is shown in Figure12 on Day PM6–PM9, which shows the distribution of total waiting time based on 10,000 simulations. According to additional calculations, the similar shape of distribution is obtained on other days and other time ranges. This variation of the total waiting time stems from traffic uncertainty effects. On this day, the estimated total waiting time based on the actual data is 316.5 minutes, but the distribution of the total waiting time between about 150 minutes and 600 minutes is observed. Even between 2.5 percentile and 97.5 percentile, the total waiting time varies between 211 and 472 minutes. This infers that the congestion level changes more than twice even if the scenario involves no uncertainty. In addition, the distribution of the total waiting time is not symmetric, and the large waiting time is observed slightly more often.Next, Figure15 shows the daily fluctuation of the total waiting time at PM6–PM9. The bar indicates the estimated total waiting time based on actual data (the same data shown in Figure 11), and the error bar indicates the 95% range of the total waiting time in the simulation. The red line shows the number of departure aircraft within this time range. On Day 15, visibility was low. From this figure, there are many interesting points found.Figure 15
Total estimated waiting time and 95% range of total simulated waiting time at PM6–PM9 on each day.First, the number of departure aircraft is almost constant in this time range, and only a small difference among the days is observed. However, like the result obtained in Table3, there is no clear relationship between the number of departure aircraft and the total waiting time. Second, even though the flight schedule is almost the same on each day in each time range, the range of the total waiting time in the simulation varies. The simulation scenario includes the initial condition of the off-block time, so it does not include the uncertainty between actual off-block time and scheduled off-block time, that is, scenario uncertainty. On Days 15 and 17, large estimated waiting time is observed, but the simulated total waiting time is also large. The low visibility seems to be a reason for large total waiting time on Day 15, but that is not the case on Day 17, maybe due to the traffic concentration. Therefore, on these days, the large total waiting time is caused not because of the traffic uncertainty effect, but because of either visibility or the scenario uncertainty, that is, the uncertainty of the off-block time compared to the scheduled time. If so, the total waiting time in the simulation varies very much with only traffic uncertainty, which means that the further large distribution will be obtained if the scenario uncertainty is considered. Therefore, the uncertainty affects the airport congestion very much, and it is almost impossible to estimate the airport congestion level on a specific day in advance, though the average congestion level can be estimated.However, this result does not directly mean that congestion cannot be relieved. If the scheduled departure or arrival time is optimally decided, the average congestion level might be reduced with keeping the traffic volume. The better spot allocation might also decrease the congestion level. However, since we understand that the uncertainty in airport traffic is significant, the uncertainty effect cannot be ignored when considering the airport operation. Otherwise, even if the taxi-out time is reduced, the uncertainty might cause the delay of take-off time or reduce the runway capacity.


## 4. Conclusions

This paper evaluated the uncertainty effect in the airport operation. The uncertainty level was obtained in each phase based on the surface movement data, and a fast-time airport traffic simulation model was developed. The validation of the simulation model was also done, and the simulation model seemed to model the uncertainty effect appropriately. Based on the developed simulation model, the characteristics of the airport traffic were investigated. The results inferred that the airport traffic congestion seemed to be mostly affected by the uncertainty of taxi-out or taxi-in and the traffic concentration, not by the traffic volume and weather conditions. Since the uncertainty effect was significant in airport traffic, it was difficult to estimate the congestion level on a specific day in advance. The developed airport simulation model would help to evaluate the airport operation with existing uncertainties.



---

*Source: 919736-2015-06-23.xml*