Skip to main content
Skip to article control options
Open AccessFull-Length Papers

Bayesian Network for Managing Runway Overruns in Aviation Safety

Published Online:https://doi.org/10.2514/1.I010726

Abstract

Runway excursions at landing constitute a major threat to aviation safety. Among them, runway overruns, defined as those occurrences when an aircraft departs the end of a runway, are the most frequent events. Although their occurrence rate is low, the entailed consequences may be very severe in terms of lives and aircraft damage. The main contributing factors to this event and their relationships are studied with the aid of a Bayesian network, while also modeling the nodes’ conditional distributions. Then, inferences and predictions are made for the quantities of interest. The issues uncovered suggest several operational recommendations to reduce the probability of facing a runway overrun when landing.

I. Introduction

According to Ref. [1], around 49% of fatal accidents involving commercial jet aircraft occur during the final approach and landing phase. One of the potential events at this stage is a runway excursion, which constitutes a major safety issue for commercial aviation systems [2]: they account for 30% of all incidents and accidents in air transport and 97% of all runway-related accidents [3]. We distinguish between occurrences in which an aircraft departs the side (veeroff) or the end (overrun) of a runway. We shall focus our attention on the latter because they are not only the most frequent events among runway excursions but they also entail the most severe consequences [4]. Reference [5] described how numerous runway lengths were designed many years ago for propeller aircraft: the growth of surrounding buildings, the advent of the turbojet aircraft, and the large increase in air traffic were not anticipated, with all of these contributing to an escalation in overruns and their impact. For these reasons, regulatory authorities and the aviation industry keep on investigating runway safety solutions. However, commercial aviation operates at such a high safety level that breakthrough improvements would require deeper analysis and more sophisticated tools and methods [68]. Runway overrun (RO) accidents occur with low probability. Yet, their implications may be severe in terms of both lives and financial costs. Thus, it is critical to further reduce their occurrence rate. Because of the potentially catastrophic consequences of ROs, novel modeling perspectives may provide new insights to enhance safety; see Refs. [911] for views on related events.

We present here a risk assessment model for ROs based on a Bayesian network (BN) [1214] aimed at helping the aviation safety community improve management of such types of threats. BNs have recently been used to deal with different safety-related aviation problems, like midair collisions or controlled flight into terrain accidents [15,16]. The goal of our paper is twofold. On the one hand, we provide a graphical model interrelating the factors identified as affecting RO with a recognized precursor of it while also describing how relevant quantities for RO risk management may be assessed. On the other hand, we populate the predictive distributions at the BN nodes based on a specific case study and make several policy recommendations.

Our ultimate intention is to provide a template that can be used as a starting point for similar applications aimed at assessing the risk of ROs in different scenarios. With the proposed model, the involved stakeholders convert their safety-related data into valuable information, providing suggestions to improve safety, as we shall show. Section II describes the features of RO events. In Sec. II, we provide our BN model for RO risk assessment. We populate the model in Sec. IV with the aid of a case study. We end up with some discussion.

II. Runway Overruns at Landing

According to Ref. [3], there were 1429 accidents worldwide entailing major (57%) or substantial (43%) aircraft damage in commercial air transport over the period of 1995–2008. Of those, almost 30% (417) were ROs: 34 of which involved fatalities. The likelihood of fatalities in such types of accidents is smaller than in other runway-related events, like incursions. However, due to their higher occurrence rate, ROs cause most aviation casualties during the landing phase. As an example, in the aforementioned period, 712 people were killed as a result of an RO. Moreover, the annual economic impact of ROs to, e.g., U.S. airlines was estimated at around 900 million U.S. dollars in Ref. [17].

In 2006, the Flight Safety Foundation (FSF) conducted a major descriptive study on runway safety [18]. After a comprehensive exploratory data analysis, they confirmed that ROs pose greater threats than other types of runway-related accidents. To mitigate such accidents, the FSF developed the “Safe Landing Guidelines” [18] summarized in Table 1.

The data suggested that the likelihood of an approach-and-landing accident increases greatly if any of the seven guidelines is not met and, even more, if two of them are not satisfied. Except for guideline d, the guidelines refer to human reliability issues (related to runway conditions at landing). The FSF recognizes guideline g as a precursor of having an RO: the risk of an accident is significantly greater if the aircraft is not slowed down to less than 80 kt by the time it reaches the point on the runway when only 2000 ft of landing distance available (LDA) remain. Typically, at 80 kt, the thrust-reverse levers are turned to reverse idle before stowing them at taxi speed.

However, to our knowledge, explicit probabilistic relationships among these variables have not yet been established. In Ref. [19], RO accident data were collected and normalized, pointing out multiple contributing factors to ROs. Reference [20] indicated that the most serious RO events could have been prevented if the Safe Landing Guidelines in Table 1 had been observed. Our approach will incorporate flight data and structured expert judgment to build a BN relating the factors potentially affecting RO events. As a consequence, we propose the adoption of several operational recommendations aimed at reducing the likelihood of having an RO.

III. Model for Runway Overrun Risk Assessment at Landing

We provide a model to facilitate the evaluation of the probability of an overrun at a runway. We intend it to be a template that can serve as a starting point for similar applications aimed to assess the risk of ROs in different scenarios. Our final goal is to compare several runways or airports in terms of their probabilities of RO at landing as a pointer to identifying more dangerous infrastructures. This may suggest airport operators, airliners, or aviation safety agencies about runways in which special care is needed under identified operational conditions. Our model will relate overrun probabilities with potentially triggering factors, suggesting the introduction of mitigation measures. We focus on assessing the probability of RO, without paying attention to the entailed consequences. Thus, we implicitly assume that we are comparing runways with a similar operational environment, as in the case in Sec. IV. Our final discussion includes pointers on consequence assessment.

We build a BN to predict ROs based on several sources: 1) the contributing factors presented in Table 1, which summarize information available from manufacturers, operators, and aviation safety agencies [18]; 2) some of the results presented in Ref. [21] related to normal operational landing performance; 3) expert judgment concerning relations between the involved variables; and finally, 4) landing data. Our initial training dataset covers 897 operations for the same airline and aircraft type (short- and medium-haul fleets); landings took place at five runways in northern Spain airports with similar operational conditions, and sometimes with tough (rain, fog, etc.) but not extreme (in terms of ice, snow, hail, etc.) weather. Their locations are not disclosed for confidentiality reasons.

As a target variable Y, we have chosen the LDR (in feet) when the aircraft speed is 80 kt, as considered in guideline g. Ignoring the presence of other variables in this discussion for the moment, we are interested in estimating Pr(Y<2000), which is used to compute

Pr(RO|Y<2000)=Pr(Y<2000|RO)Pr(RO)/Pr(Y<2000)

The numerator probabilities may be estimated from data in the FSF studies mentioned previously based on standard beta-binomial models [22]. The other nine variables included in our BN are listed in Table 2, and they are related to those appearing in the FSF guidelines [18].

We used GeNIe’s [23] built-in Bayesian search learning (BSL) algorithm [24] to build our BN. A first version of it was provided by our experts, which served as the initial solution to the algorithm. Because GeNIe does not support continuous variables, we discretized them using its built-in hierarchical algorithm [25]. We then fed BSL to obtain our solution using GeNIe’s default clustering algorithm for inference [26]. Additional information about learning BN structures can be found in Refs. [27,28].

Several experts checked the soundness of the suggested variable relations in the network in Fig. 1.§ The experts found all (appearing or missing) links pertinent, except for the possible connection between cross and tailwind, which did not seem relevant according to the used data, as confirmed by the slightly higher Bayesian information criterion (BIC) [29] network score computed when removing that link. We decided to leave it as such and check for possible correlations in specific cases, as in Sec. IV. Then, the proposed BN defines a probabilistic model relating factors connected with ROs through

p(a,b,c,d,e,f,g,h,i,y)=p(a)p(b|a)p(c|a)p(d|a)p(e|a,b,c)p(f|a,e)p(g|a)p(h|a,e)p(i|e)p(y|a,d,f,g,h,i)(1)
where we use p to designate a generic distribution.

Fig. 1
Fig. 1

BN of factors related to excursion events (Rwy, runway; Xwind, crosswind; Rwy_cond, runway condition; Unst_app, unstabilized approach; T_rev_full, time reverse full; Height_thr, height threshold; Diff_IAS_Vapp_thr, Difference between IAS and Vapp at threshold; Rmg_rwy_80, remaining runway at 80 kt).

We now discuss several potential model uses derived from Eq. (1), which are later illustrated in Sec. IV. As mentioned, for a given organization, we shall compute the target probability Pr(Y<2000) to support the assessment of its probability of RO. We may use it to rank airliners, airports, or runways under similar environmental (weather at landing or state of the runway) and operational conditions. This probability may be obtained through

Pr(Y<2000)=02000p(y)dy=02000[aΩadΩdΩfgΩgΩhΩieΩeΩbΩcp(y|a,d,f,g,h,i)×p(a)p(d|a)p(f|a,e)p(g|a)p(h|a,e)p(i|e)p(e|a,b,c)p(b|a)p(c|a)dcdbdidhdf]dy(2)
which is typically approximated through Monte Carlo (MC) simulation with
1N#{k:1kN,yk<2000}
for a sample
{yk}k=1Np(y)
where the # symbol means the “number of”; and Ωa,,Ωi denote the domains of variables A to I, respectively. We may use the composition sampler in Algorithm 1 to obtain the sample:

The target probability at a given runway a̲ is obtained using Bayes’s formula

Pr(Y<2000|A=a̲)=Pr(Y<2000,A=a̲)/Pr(A=a̲)
The numerator would be approximated through
1N#{k:1kN,yk<2000}
where
{yk}k=1N
is a sample from Algorithm 1 fixing ak=a_. The factor Pr(A=a̲) would be obtained from the initial assessments.

We are also interested in computing the aforementioned probabilities under certain operational conditions, as we will outline in the following. Such quantities could provide crews with crucial warnings about potentially dangerous situations. For brevity, we omit most computational details because they resemble those provided for Pr(Y<2000|A=a̲), and we relegate the technical aspects of the following discussion to Appendix A.

We first assess the target probability at a given runway a_ under strong wind conditions [say, Pr(Y<2000|A=a̲,B>b̲,C>c̲)] for threshold tail and crosswind speeds b_ and c_. This helps in setting alarms to the crew, letting them know an eventual important increase in RO probability, given the wind conditions forecast. A similar approach may be used to assess the influence of the reverse thrust (RT) in runway a̲ by computing Pr(Y<2000|A=a̲,Ff̲) for a certain cutoff maximum reverse thrust (MRT) time f_.

Another key aspect for operators is the influence of stability, which is assessed through the probability

Pr(Y<2000|A=a̲,E=Unst)/Pr(Y<2000|A=a̲,E=Stab)

It has a likelihood ratio form comparing risks at a given runway a̲, depending on the approach stability. Should this be bigger than one, we would alert about dangerous conditions associated with an unstabilized approach, suggesting the implementation of procedures to increase situational awareness.

It is also relevant to determine the influence of runway conditions through

Pr(Y<2000|A=a̲,D=Contaminated)/Pr(Y<2000|A=a̲,D=Noncontaminated)

When the ratio is much bigger than one, it suggests the strong influence of runway conditions, pointing out toward improving situational awareness among crew if the runway is contaminated.

IV. Template for Case Studies

We present a case study that may be used as a template for related applications aimed at assessing the risk of ROs in different scenarios with similar operational and environmental conditions. It is undertaken by an airline willing to assess RO risk at three runways designated rwy-01, rwy-02, and rwy-03. Within the framework of RO risk management, the proposed graphical model (Fig. 1) and the models and priors at various nodes identified in the following could be used as starting points for other case studies considered, e.g., for other airports, runways, companies, or countries. The case study also serves to show the operational insights and recommendations that may be derived from our model.

We have available a second dataset provided by an operator, which is independent of that used in Sec. III to learn the BN structure. We use it, together with expert judgment, to build the predictive distributions at the network nodes. References [30,31] provided examples in which solely expert judgment was used in building and populating a BN model for, respectively, highway and nuclear plant safety assessments. In turn, Refs. [32,33] illustrated how to learn the BN structure and parameters solely from data. Due to the nature of the RO problem, we found very few critical cases in the period considered. However, we used expert opinion as another source of information to provide forecasts and combined it with the little data available. Reference [34] included additional information about learning the node distributions. Some of the priors we use are diffuse but proper (i.e., with large variance but integrating to one); others reuse knowledge gained from data in other parts of the problem; finally, some other priors model the experts’ knowledge in relation with the physical aspects of the problem. In any case, they are object to a sensitivity analysis in Sec. IV.B.

We therefore used data from n=266 landings over a period of 10 months, with all referring to the same airline and aircraft type. They concern the three runways mentioned previously and shown in Fig. 2 (all with LDAs smaller than 7200 ft), as well as similar operational and environmental conditions around the airports. In these runways, should an RO happen, it would just affect the aircraft and/or passengers and crew because there are no buildings nearby. Thus, we pay no attention to consequences because they would be similar for the three involved infrastructures (in terms of fatalities, injuries, and aircraft damage or destruction), and we focus on estimating the probability of the relevant event. See our final Sec. V for a discussion on consequence assessment.

Recalling our analysis in Fig. 1, we start by investigating possible relationships between cross- and tailwind components in this case. Pearson’s correlation coefficients are, respectively, r1=0.12, r2=0.04, and r3=0.21 if we include zero tailwinds; after removing such null values, the correlations become r1=0.04, r2=0.01, and r3=0.35. None suggest strong dependence. Therefore, we proceed with our graphical model from Fig. 1 with no direct link between cross- and tailwind nodes.

We illustrate now several model uses supporting relevant operational decisions. We describe in Appendix B the predictive distributions at various nodes deployed to populate model (1). For simplicity, we drop the dependence on data, writing, e.g., p(A) rather than p(A|data). In most nodes, we shall be able to obtain only a sample from the relevant predictive distribution. When pertinent, we display the expected predictive distribution over the data histogram to suggest the appropriate model fit. All computations were performed with R and WinBUGS [35,36].

A. Using the Model

We have developed a sampling scheme to implement the computations in Sec. III. The reported results were obtained after 1) generating samples of a size of 1,000,000 from the involved predictive distributions of the variables A,B,,I,Y, which are detailed in Appendix B; 2) performing the relevant calculations (probabilities and likelihood ratios) as outlined in Sec. III; and 3) replicating such a scheme 100 times to estimate the corresponding uncertainty.

Fig. 2
Fig. 2

Runways for the case study (source: Google Earth Pro).

First, we are interested in the overall probability of having a landing in which the landing distance remaining (LDR) is less than 2000 ft when the aircraft is at 80 kt. Through Algorithm 1, we have obtained Pr^(Y<2000|data)=0.00019 with an associated standard deviation of 4.43105. Thus, given the load of the incumbent operator, we should expect one landing with high risk of RO at any of the three runways approximately every 400 days. Extrapolating this result to the operations of the airline in its country where it operates (assuming all runways have similar operational environments), we should expect one hazardous operation (in relation to RO) approximately every 26 days. If this is deemed excessive by the airline, this would require risk mitigation strategies.

We find out which of the three runways is the most dangerous one. Using Algorithm 1 with ak=a̲, we get the estimates in Table 3.

Rwy-01 has the largest probability, entailing roughly one hazardous operation every two years. Figures are comparable for rwy-02, whereas the target probability at rwy-03 is approximately 10 times smaller. Rwy-01 is indeed the most complicated runway from an operational viewpoint because 1) it is not equipped with an instrument landing system (ILS) and pilots have to approach it visually; 2) it has a slightly U-shaped valley form; and 3) adverse wind conditions are more frequent (see Figs. 3 and 4). Eventually, risk management in such a runway would proceed by introducing an ILS and emphasizing wind issues in briefing manuals and arrival charts.

Fig. 3
Fig. 3

Crosswind component at threshold by runway, with overlaid predictive densities.

Fig. 4
Fig. 4

Positive tailwind factored by runway, with overlaid predictive densities.

In regard to the influence of wind conditions, we consider b̲{10,15,20} and c̲{0,5,10} as relevant cutoff wind speeds and assess the target probabilities in Table 4. We use the samples generated in Algorithm 1 with ak=a̲ for such a purpose.

Clearly, higher probabilities hold as wind conditions worsen: landing under such circumstances tends to increase the probability of an unstabilized approach. This will affect the MRT, the difference I, and the height at threshold (see Fig. 1), eventually increasing the risk of having less LDR at 80 kt. If that is the case, the autobrake system must be selected. Nevertheless, under severe wind conditions (e.g., b̲=20  ktandc̲=10  kt), “going around” and proceeding to the alternate airport would be the right decision to avoid putting safety in jeopardy. This should be briefed to crews.

We analyze now the influence of RT. Its use is linked with policies to reduce fuel consumption and noise at airports: pilots are requested to operate RT at idle during ground roll, except when safety reasons recommend it. How long pilots operate at MRT will depend on the speed the aircraft needs to be slowed down, which is influenced by factors such as tailwind. RT cutoff values were set at f̲{0,5,10}. Using the samples from Algorithm 1 with ak=a̲, we obtain Table 5.

Therefore, the longer the MRT has been operated, the higher the target probability. This indicates that, if RT has been used for too long (say, over f̲=10  s), the aircraft might have experienced difficulties in decelerating during ground roll, risking an RO. Thus, this variable could be used as an RO risk precursor: relevant information for flight data monitoring teams.

We consider now the influence of approach unstability. After running Algorithm 1 with ak=a̲ and ek=e̲, we have obtained ratios (9.11, 11.56, 2.59) for the runways, respectively. Thus, more dangerous conditions will arise under an unstabilized approach, especially at rwy-01 and rwy-02. Recommendations oriented toward increasing situational awareness should be given to crews. Innovative mitigation measures could be a notification from the flight dispatcher to the crew of such high risk level or including a warning in briefing manuals or arrival charts.

Finally, we discuss the influence of runway conditions. Through Algorithm 1 with ak=a̲ and dk=d̲, we obtained ratios (1.08, 0.96, 0.95) for the three runways. They are relatively close to one, suggesting little influence of surface contamination on the risk of runway overrun in this case: possibly because pilots usually take extra precautions when facing adverse runway conditions. Thus, no specific information needs to be included in briefing manuals or arrival charts.

B. Sensitivity Analysis

We have performed a sensitivity analysis over various priors to check for robustness of the results. We have focused on those parameters with larger posterior variances to ascertain the influence of different prior choices. As an illustration, and because of its judgmental nature, we explore here sensitivity to prior variations in relation with height at the threshold. We tried alternative priors for the w, ν, and μ, reflecting different precision regarding prior knowledge. We considered the eight combinations of [ϕ,θ,(α,β)] in Table 6.

The results obtained did not show significant changes, suggesting robustness: the estimates for the corresponding predictive distribution varied by less than 5%, whereas the target probabilities differed by less than 2% with respect to those in Table 3. Similar analyses were performed for other model parameters, confirming the robustness of our priors and results.

V. Conclusions

Although they occur with very low probability, runway overruns constitute a major threat to aviation safety because their implications may be very severe in terms of lives and destroyed aircraft. Several international aviation organizations like the FSF, the International Civil Aviation Organization, or aircraft manufacturers are launching efforts to reduce landing overruns, investigating various safety strategies.

As cited, information available from manufacturers, operators, and safety agencies (The Boeing Company, the FSF, the Federal Aviation Administration, Honeywell, Airbus, etc.); expert judgment; and data were used to propose a BN defining a probabilistic model for assessing RO probabilities. The critical precursor event considered was whether the LDR at 80 kt was less than 2000 ft. Its probability was used to assess the RO risk and rank various aviation systems like runways, airports, and others. The methodology is general, and the probabilistic graph underlying the current network may be used in other cases, except for the possible relation between the crosswind and tailwind, which would need to be checked as explained. The node predictive distributions would also need to be reevaluated in other case studies, but the current proposals may serve as initial models when comparing homogeneous aviation systems. The current working BN can easily integrate new data, accordingly updating the node predictive distributions; see Ref. [34] for a discussion.

As with many other statistical models, one of the main limitations of the current proposal might be the scarcity of data. However, in this case, this is a fortunate circumstance because it means that few accidents have actually occurred. Nevertheless, it is believed that it is not a matter of having enough data but of using them properly in combination with the expert judgment available. Furthermore, as data accumulate, the corresponding predictive distributions would be updated.

The validity of the current methods has been thoroughly assessed. First, the BN built in Fig. 1 was chosen among a set of suitable models as the one with the lowest BIC. Second, each of the predictive distributions was validated for the goodness of fit, displaying (when pertinent) the expected predictive distribution over the data histogram. Some assessments were also made based on the predictive probabilities obtained. To wit, the results in Table 3 were in accordance with the observed data where, for instance, only one out of 80 operations took place under an unstabilized approach at rwy-02.

From the current case study, several general relevant issues were uncovered, including the following:

  1. The first relevant issue involves the weather conditions. Landing under windy conditions tends to increase the probability of an unstabilized approach.
  2. The second relevant issue involves the variables with the strongest influence. The variables having stronger influence over the LDR at 80 kt (and hence over RO) were the autobrake system and the difference I at threshold.
  3. The third relevant issue involves ducking under. The presence of a ducking under effect was alerted against, which is not always a good procedure.
  4. The fourth relevant issue involves unstabilized approach. As RO accident reports highlight, much more danger arises from unstabilized approaches.
  5. The fifth relevant issue involves reverse thrust. Excessive use of MRT suggests difficulties for decelerating during ground roll.

Note that the focus has been on evaluating the probability of a leading precursor of RO because homogeneous runways are being dealt with. Thus, attention has not been paid to forecasting the entailed consequences, which would typically include the number of deaths and injured people, the potential destruction of the aircraft, the induced delay, and so on. Forecasting them would require additional modeling efforts. They would then be evaluated through a multiattribute loss function; see Refs. [37,38] for the details of a general framework for aviation safety risk management. With such ingredients, the acceptability of certain risks could be assessed or appropriate countermeasures (possibly referring to signals, markings, lighting, or updated procedures that could potentially reduce overruns) could be identified, minimizing the expected loss associated with their introduction.

Many safety-related issues are yet to be explored in the aviation industry. To cite just a few, the potential impact of human behavior on the performance when landing, the influences of other external factors (such as, e.g., the congestion of airports or the pressure of tight schedules), or the attrition or lack of maintenance of some aircraft are still, to a great extent, open questions. New methodological solutions, possibly following similar approaches to the one in this paper, would be necessary to gain better understanding and propose innovative solutions to such problems.

M. J. KochenderferAssociate Editor

§ GeNIe uses the thickness of the arcs to indicate the strength of influence between the nodes they connect. For ease of viewing, we have further emphasized stronger intensities with darker links.

¶ In what follows, we will emphasize in italics the risk management options suggested.

Appendix A: Probabilities and Likelihood Ratios

A.1. Computation of Pr(Y<2000|A=a̲,B>b̲,C>c̲)

Using Bayes’s formula, and canceling common terms, this probability is proportional to

1N#{k:1kN,yk<2000,bk>b_,ck>c_}
based on Algorithm 1 with ak=a_, using the output {yk,bk,ck}k=1N. The proportionality factor arising from Bayes’s formula, Pr(B>b̲|a̲)Pr(C>c̲|a̲), is obtained from the airport assessments on the crosswind and tailwind.

A.2. Computation of Pr(Y<2000|A=a̲,E=Unst)/Pr(Y<2000|A=a̲,E=Stab)

After canceling common terms, the numerator Pr(Y<2000,A=a̲,E=e̲) is approximated through

1N#{k:1kN,yk<2000}
where {yk}k=1N is obtained from Algorithm 1, fixing ak=a_ and ek=e_. The denominator is approximated through
Pr(E=e̲|A=a̲)1Nk=1NPr(E=e̲|a̲,bk,ck)
for a sample {bk,ck}k=1N from Algorithm 1 with ak=a_.

A.3. Computation of Pr(Y<2000|A=a̲,D=Contaminated)/Pr(Y<2000|A=a̲,D=Noncontaminated)

The required probability is proportional to

1N#{k:1kN,yk<2000}
with {yk}k=1N as a sample from Algorithm 1 fixing ak=a_ and dk=d_. The proportionality factor Pr(D=d̲|a̲) is obtained from the initial assessments.

Appendix B: Predictive Models at Nodes

Runway usage p(A) is a categorical variable with values A{1,2,3} referring to the runway used. We employ a Dirichlet-multinomial model with a uniform prior for runway usage [22]. Let (p1a,p2a,p3a) be the respective usage proportions of such runways. Then,

(p1a,p2a,p3a)Dir(1,1,1),(x1a,x2a,x3a)M(n;p1a,p2a,p3a),(p1a,p2a,p3a)|dataDir(1+x1a,1+x2a,1+x3a)(B1)
with x1a=104, x2a=80, and x3a=82 as the number of operations at each runway. When needed, we use the predictive means
Pr^(A=|data)=(1+xa)/(3+n),=1,2,3
as estimates of runway usage. In this case,
Pr^(A={1,2,3}|data)=(0.39,0.30,0.31)
suggesting similar usage for all runways.

We model now the crosswind, given runway p(B|A). All 266 flights entailed nonzero crosswind, with histograms in Fig. 3. For each runway, we use a gamma model B|A=G(νb,νb/μb) [39].

We respectively consider diffuse, but proper, exponential νbE(0.01) and inverse gamma μbIG(1,1) priors. There is no closed-form expression for the posteriors, but we may sample from μb|data and νb|data, as summarized in Table B1. Based on such samples, we approximate the predictive distribution p(b|A=,data), plotted with a dashed line in Fig. 3, suggesting a good fit.

Next, we model the tailwind, given runway p(C|A). Out of the 266 flights, 125 took place with no tailwind (Table B2). Figure 4 presents the histograms for the cases with positive tailwind.

For each runway , we model the tailwind distribution C through a mixture

C=p0cI0c+p1cC+,p0c+p1c=1,p0c,p1c0(B2)

I0c describes zero tailwind cases, and C+ describes the positive ones. Based on a uniform prior for the beta-binomial model [22], we get that, a posteriori,

p1c|dataBe(1+xc,1+nxc)(B3)
where n is the number of operations at runway , and xc is the number of operations with positive tailwind; we may estimate p1c with its posterior mean (1+xc)/(2+n). We model the positive tailwind component through a gamma model, as for the crosswind C+G(νc,νc/μc), with diffuse but proper priors. Table B3 summarizes the posteriors. The predictive distribution p(c|A=,data) is approximated through MC simulation and displayed in Fig. 4, suggesting again a good fit.

We deal now with modeling p(D|A), referring to runway conditions at A, given the runway. We have available data about weather conditions for the incumbent operations (0 denotes noncontaminated, and 1 denotes contaminated). We consider a runway contaminated if rainfall over the previous 24 h exceeded 8  mm/m2 [40]. We use a beta-binomial model with a uniform prior for p1d, with the probability of runway being contaminated when landing. A posteriori, p1d|dataBe(1+xd,1+nxd), where xd is the number of operations on a contaminated runway. When necessary, we summarize p1d through its posterior mean (1+xd)/(2+n). Table B4 shows the estimates and the available data in parentheses. The predictive distribution is [41]

p(d|A=,data)BetaBin(1,1+xd,1+nxd)

Then, we model stability given the runway, crosswind, and tailwind through p(E|A,B,C). We label ej=1 (0) if the approach is unstabilized (stabilized). We use a logistic regression model for the probabilities pje=Pr(Ej=1|A=,B,C):

logit(pje)=log(pje/(1pje))=αe+β1eBj+β2eCj,j=1,,xa,=1,2,3(B4)
where Bj and Cj are, respectively, the jth values of the crosswind and tailwind components at threshold when landing at runway . We assume diffuse but proper normal priors N(0,105) for the α, and we assume E(1) for β, because we expect greater probabilities of unstabilized approaches under worse wind conditions. We obtained the posterior means and standard deviations in Table B5. Note that the posterior estimates of β are slightly positive in all runways, although the crosswind β1e tends to have greater influence over approach stability than tailwind β2e. As suggested by our experts, this could be due to the fact that crews would be more aware when landing under a strong tailwind, actually mitigating its effect. However, the corresponding posterior variances are relatively large, especially that of β2e at rwy-02 (for which only one operation under the unstabilized approach was available), suggesting uncertainty given the little data available.

Because there is no closed-form expression for the predictive distribution, we use a composition simulation scheme to obtain the required predictive samples. Symmetric 0.9 predictive intervals [22] covered observations, suggesting good performance of our model.

We focus now on modeling RT, given the runway and stability p(F|A,E). Out of 266 operations, 118 entailed reverse at idle thrust. The remaining 148 landings used MRT during ground roll. Table B6 summarizes these data, factored by landing runway and approach stability. Histograms are shown in Fig. B1.

Fig. B1
Fig. B1

MRT factored by approach and runway, with overlaid predictive densities.

We distinguish between stabilized and unstabilized approaches. For stabilized ones, we decompose RT F through the mixture

F=p0fI0f+p1fF+,p0f+p1f=1,p0f,p1f0(B5)
for each runway: =1,2,3. The first component I0f describes the idle RT cases, whereas F+ models the MRT ones. We use the same modeling approach as for tailwind. Posteriors are summarized through means and standard deviations in Table B7.

For unstabilized cases, there are no zero values (Fig. B1). We model the MRT through a gamma model, using as priors the posterior distributions obtained for the positive cases in the stabilized approach, i.e.,

F|(A=,E=1)G(νf,u,νf,u/μf,u),νf,uνf,s|data,μf,uμf,s|data(B6)
where =1,2,3 represents the runway. Posteriors are summarized in Table B8.

We proceed as with the tailwind to obtain a sample from the predictive densities. Overlaid expected predictive densities in Fig. B1 suggest a good fit.

We consider autobrake given the runway p(G|A). We use a Dirichlet-multinomial model with a uniform prior for each runway. Based on the coding (1 denotes no, 2 denotes low, and 3 denotes medium) autobrake, we have

(p1g,p2g,p3g)Dir(1,1,1),(x1g,x2g,x3g)M(n;p1g,p2g,p3g)(B7)
where (p1g,p2g,p3g) represent the usage proportion of the three autobrake system configurations at runway , and (x1g,x2g,x3g) are the corresponding numbers of operations under such configuration. A posteriori,
(p1g,p2g,p3g)|dataDir(1+x1g,1+x2g,1+x3g)(B8)

When needed, we use estimates

Pr^(G=γ|A=,data)=(1+xγg)/(3+n),γ=1,2,3

See Table B9. Note that the autobrake system is mostly used in its “medium” setting at the shortest runway (rwy-01) because it provides the highest deceleration.

Next, we model p(H|A,E), referring to height at the threshold, for a given runway and stability conditions. Figure B2 represents the height H when aircraft overfly the threshold. Ideally, this should be done at 50 ft (vertical dashed line). For each runway, we consider that data come from a mixture of, at most, three gamma distributions: one for values close to zero (too low), and the other two for intermediate to upper values (too high),

H|(A=,E=0)m=13wmh,sG(νmh,s,νmh,s/μmh,s),wmh,s0,m=13wmh,s=1(B9)

Fig. B2
Fig. B2

Height at threshold factored by approach and runway, with overlaid predictive densities.

We consider the prior

(w1h,s,w2h,s,w3h,s)Dir(ϕ),νmh,sE(θ),μmh,sIG(αm,βm)(B10)
for m=1,2,3, with a restriction of μ1h,s<μ2h,s<μ3h,s to ensure identifiability. Because we have little information about the mixture weights, we set ϕ=(1,1,1) to provide a uniform prior over them. Matching our expert’s prior knowledge, we reflect that the first group of heights will have a mean of around 10 ft (for pilots crossing the threshold too low); the second around 40 ft (slightly below the recommended 50 ft at threshold); and, finally, the third one around 60 ft (slightly over 50 ft), choosing (α1,α2,α3)=(27,66,38) and (β1,β2,β3)=(260,2600,2220) as the parameters of the inverse gamma prior distributions for μm, with expected values of (10, 40, 60) and standard deviations of (2, 5, 10) accounting for the expert’s greater uncertainty at higher altitudes. Finally, we choose a relatively diffuse (i.e., with relatively large variance) exponential prior for νm, m=1,2,3 by setting θ=0.05. As there is no closed-form expression for the posteriors, we use a Markov Chain Monte Carlo simulation. The obtained estimates are in Table B10. We cannot explicitly obtain the predictive distribution p(h|A=,E=0,data), which we approximate with a sampling scheme.

For unstabilized approaches, pilots tend to fly either slightly below or above 50 ft, according to experts. Therefore, we consider only two mixture terms:

H|(A=,E=1)m=12wmh,uG(νmh,u,νmh,u/μmh,u)(B11)

Due to the small sample size, we use as priors the posteriors over the last two components in the stabilized case and their corresponding (renormalized) weights. Posterior estimates are shown in Table B11.

For the predictive p(h|A=,E=1,data), we proceed through a sampling scheme. The overlaid expected predictive densities in Fig. B2 suggest a good fit.

We deal now with modeling IASVapp, given the stability of approach p(I|E). Figure B3 represents I=IASVapp at threshold, which is factored by the type of approach E. Ideally, I=0; i.e., the indicated airspeed (IAS) should coincide with Vapp, which is indicated with a vertical dashed line.

Fig. B3
Fig. B3

Difference IASVapp factored by approach, with overlaid predictive densities.

For both types of approaches, we assume a normal model I|EN(μ,σ2) with a noninformative prior p(μ,σ2)1/σ2. A posteriori,

p(μ,σ2|data)σne2exp(12σ2[(ne1)se2+ne(μi¯e)2])(B12)
where ne is the number of operations; i¯e is the sample average of (IASVapp); and se2 is the sample variance, with e=0 (1) indicating an (un)stabilized approach. Note that p(μ|data) is a Student’s t distribution with (ne1) df, location i¯e, and scale parameter se2/ne [42]. The predictive p(i|E=e,data), which is another Student’s t distribution with the same df and location as p(μ|data) but a different scale parameter (1+1/ne)se2, is overlaid in Fig. B3, suggesting a good fit. Table B12 summarizes the parameters.

Finally, we model the LDR Y at 80 kt, given the relevant predecessors p(Y|A,D,F,G,H,I). We plot in Fig. B4 the LDR when the aircraft is decelerated to 80 kt for each runway, indicating with a vertical dashed line the threshold distance of 2000 ft.

Fig. B4
Fig. B4

Remaining runway factored by runway, with overlaid predictive densities.

We use a linear regression model

yj=αy+β1yDj+β2yFj+β3y,lGjl+β3y,mGjm+β4yHj+β5yIj+ϵj(B13)
for j=1,,xa, =1,2,3, where Fj, Hj, and Ij are, respectively, the jth value of MRT, the height at threshold, and the difference between the IAS and Vapp when landing at runway ; Dj describes the runway condition D; Gjl and Gjm refer to the autobrake variable G: gjl (gjm) is one if the autobrake is set to low (medium), and it is zero otherwise; finally, ϵj is the error term, which we assume is Student’s t distributed to achieve greater robustness, with a zero location vector and a scalar dispersion matrix σϵ2 [43]. We assume a diffuse but proper joint prior for σϵ2 and the vector β of regression coefficients p(β,σϵ2|)=p(σϵ2|)p(β|σϵ2,), where σϵ2 follows a Fischer’s distribution F(1,1), and β|σϵ2 is modeled as a Student’s t distribution tν(0,τ), with ν (the degrees of freedom) following (in turn) a discrete uniform distribution in [2, 100] and τ (the precision, inverse of the variance) as an inverse gamma IG(0.01,0.01). Posterior summaries appear in Table B13.

The LDR Y at 80 kt seems strongly influenced by the autobrake G and the difference I at threshold. The positive value of β3y,m suggests that when, the autobrake is selected at medium, the landing roll will decrease, whereas the LDR will increase. The same reasoning holds for β3y,l, albeit to a lesser extent. With respect to the negative sign of β5y, this suggests that a greater IAS at threshold implies greater kinetic energy and, as a consequence, more runway is required. As for the runway condition D, β1y has large positive values for rwy-02 and rwy-03. This could be due to two factors:

  1. The runway is severely contaminated at landing (typically with heavy rain), which is a relatively frequent phenomenon in the incumbent locations, especially in winter; this might damp the airplane during ground roll, therefore reducing RO risk.
  2. When pilots observe adverse weather conditions at destination, they will usually increase their situational awareness, actually leading to safer landings.

The influence of runway conditions seems smaller for rwy-01, possibly suggesting that, in some cases, runway contamination might eventually cause hydroplaning, reducing the efficacy of the deceleration process. Nevertheless, the posterior variance of β1y implies nonnegligible uncertainty, especially for rwy-02 and rwy-03. The impact of the other variables seems less relevant. For F, because β2y>0, this suggests that the LDR at 80 kt is slightly increased by the use of MRT. For H, β4y>0 could be an indicator of the tendency of pilots to pitch down before the threshold so as to touch down as soon as possible, a maneuver called ducking under: pilots sometimes descend below the glide slope to have a larger LDA, especially at short runways. However, when ducking under is performed incorrectly, it could have an opposite effect due to the increase in kinetic energy and the ground effect, as well as the reduction in induced drag. Thus, pilots should be advised to avoid an excessive or inappropriate use of ducking under in these runways.

A closed-form expression for the required predictive distribution is not available, but we may draw from the runway predictive distribution. Symmetric 0.9 predictive intervals covered observations, suggesting good performance of our model.

Acknowledgments

This work is supported by the Spanish Ministry of Economy and Innovation project MTM2017-86875-C3-1-R and the AXA-ICMAT Chair on Adversarial Risk Analysis. Data concerning weather conditions at the case study runways have been provided by the Spanish Agencia Estatal de Meteorologa.

References

Tables

Table 1 Safe landing guidelines

GuidelineDescription
aFly a stabilized approach.
bHeight at threshold crossing should be at most 50 ft.a
cSpeed at threshold crossing to be within 10 kt of the reference airspeed Vref.
dTailwind should be no more than 10 kt (00 kt) for noncontaminated (contaminated) runways.
eTouch down on runway centerline at touchdown point.
fAfter touchdown, promptly transition to the desired deceleration setting.
gSpeed should be less than 80 kt with 2000 ft of landing distance remaining.

aRunway thresholds are markings across the runway, delimiting the beginning and end of the designated space for landing and takeoff under standard conditions.

Table 2 Variables in the BN

VariableDefinition
ARunway. Categorical: identifies the relevant runway.
BCrosswind (knots). Nonnegative continuous, discretized to nearest integer.
CTailwind (knots). Same as crosswind.
DRunway condition. Categorical: contaminated or not.a
EType of approach. Categorical: stabilized or not.
FMaximum reverse thrust (seconds). Continuous. Time of MRT operation during ground roll.b
GAutobrake. Categorical: no autobrake, low and medium.
HHeight at threshold (feett). Nonnegative continuous; equals zero when aircraft touches down.
IIASVapp (knots). Continuous. Difference between IAS and Vapp at threshold.c
YLDR at 80 kt (feet).

aRunway contamination may be due to factors such as rain, standing water, hail, ice, or snow.

bNote that the type of large commercial jet aircraft we have analyzed (short- and medium-haul fleets) is typically equipped with such a system.

cIAS is the airspeed read directly from instruments. Vapp represents the airspeed to be maintained down to 50 ft over the runway threshold plus corrections for wind.

Algorithm 1: Sampling scheme for the target variable

For k=1 To N
1:Sample akp(a)
2:Sample bkp(b|ak), ckp(c|ak), dkp(d|ak), gkp(g|ak)
3:Sample ekp(e|ak,bk,ck)
4:Sample fkp(f|ak,ek), hkp(h|ak,ek), ikp(i|ek)
5:Sample ykp(y|ak,dk,fk,gk,hk,ik)
6:k=k+1

Table 3 Posterior probabilities Pr(Y<2000|A=a̲,data) with standard deviations (in parentheses) (×104)

 Data
Rwy-013.27 (0.56)
Rwy-022.47 (0.51)
Rwy-030.29 (0.84)

Table 4 Posterior probabilities Pr^(Y<2000|A=a̲,B>b̲,C>c̲,data)(×104)

 Rwy-01Rwy-02Rwy-03
b̲=10b̲=15b̲=20b̲=10b̲=15b_=20b̲=10b̲=15b̲=20
c_=07.819.7211.081.791.942.140.210.340.49
c_=59.4310.5712.111.982.172.310.260.390.48
c_=1011.8412.8113.442.773.954.720.600.850.99

Table 5 Posterior probabilities Pr^(Y<2000|A=a̲,Ff̲,data)(×104)

f_=0f̲=5f_=10
Rwy-01
0.361.311.72
Rwy-02
0.180.540.84
Rwy-03
0.00180.00770.016

Table 6 Sensitivity analysis

wDir(ϕ)νE(θ)μIG(α,β)
ϕθ(α,β)
(0.4,0.4,0.2)0.1((10,10,10),(200,400,600))
(0.5,0.4,0.1)0.2((2.2,2.2,2.2),(40,50,100))

Table B1 Posterior summaries of parameters for B|A with standard deviations (in parentheses)

 Rwy-01Rwy-02Rwy-03
νb|data3.76 (0.50)3.96 (0.60)6.11 (0.91)
μb|data6.50 (0.34)5.21 (0.29)4.03 (0.18)

Table B2 Presence of tailwind by landing runway

TailwindRwy-01Rwy-02Rwy-03
No492947
Yes555135

Table B3 Posterior summaries of parameters for C|A with standard deviations (in parentheses)

 Rwy-01Rwy-02Rwy-03
νc|data2.29 (0.40)3.16 (0.58)2.78 (0.62)
μc|data4.49 (0.41)4.34 (0.35)4.24 (0.44)

Table B4 Posterior estimates and data for runway conditions factored by runway

Runwwy conditionsRwy-01Rwy-02Rwy-03
Noncontaminated (0)0.89 (93)0.93 (75)0.89 (74)
Contaminated (1)0.11 (11)0.07 (5)0.11 (8)

Table B5 Posterior summaries of parameters for p(E|A,B,C) with standard deviations (in parentheses)

 Rwy-01Rwy-02Rwy-03
αe|data10.96 (3.81)8.33 (3.02)8.47 (2.64)
β1e|data0.54 (0.26)0.49 (0.29)0.59 (0.33)
β2e|data0.39 (0.23)0.11 (0.37)0.48 (0.23)

Table B6 RT during ground roll

 RT
 StabilizedUnstabilized
 Rwy-01Rwy-02Rwy-03Rwy-01Rwy-02Rwy-03
Idle394137000
Maximum633842213

Table B7 Posterior summaries of parameters for F|(A,E=Stab) with standard deviations (in parentheses)

 Rwy-01Rwy-02Rwy-03
νf,s|data5.15 (0.88)6.16 (1.35)3.35 (0.70)
μf,s|data9.47 (0.52)10.24 (0.70)9.45 (0.83)

Table B8 Posterior summaries of parameters for F|(A,E=Unst) with standard deviations (in parentheses)

 Rwy-01Rwy-02Rwy-03
νf,u|data5.27 (0.88)3.95 (0.84)3.34 (0.66)
μf,u|data9.46 (0.52)9.98 (0.82)9.39 (0.79)

Table B9 Posterior estimates and data for autobrake factored by runway

AutobrakeRwy-01Rwy-02Rwy-03
No (1)0.21 (22)0.38 (31)0.51 (42)
Low (2)0.12 (12)0.28 (22)0.27 (22)
Medium (3)0.67 (70)0.34 (27)0.22 (18)

Table B10 Posterior parameters for gamma mixture (stabilized approach) with standard deviations (in parentheses)

 Component 1Component 2Component 3
Runwayν1^h,sμ1^h,sw1^h,sν2^h,sμ2^h,sw2^h,sν3^h,sμ3^h,sw3^h,s
Rwy-012.00 (0.55)9.98 (1.65)0.33 (0.07)14.03 (11.70)40.13 (10.49)0.27 (0.19)13.19 (10.65)63.37 (12.83)0.40 (0.20)
Rwy-021.74 (0.40)9.68 (1.42)0.51 (0.09)23.71 (18.54)26.96 (8.22)0.19 (0.14)27.40 (14.29)44.45 (12.37)0.30 (0.16)
Rwy-031.76 (0.48)10.54 (2.13)0.40 (0.12)9.57 (10.80)33.38 (9.04)0.32 (0.19)16.24 (20.40)68.64 (18.12)0.28 (0.19)

Table B11 Posterior parameters for gamma mixture (unstabilized approach) with standard deviations (in parentheses)

 Component 1Component 2
Runwayν1^h,uμ1^h,uw1^h,uν2^h,uμ2^h,uw2^h,u
Rwy-0114.72 (11.23)39.77 (10.21)0.67 (0.23)12.88 (10.33)63.89 (13.21)0.33 (0.21)
Rwy-0224.11 (18.21)27.12 (8.09)0.71 (0.14)27.76 (14.93)43.94 (11.99)0.29 (0.16)
Rwy-039.03 (9.89)33.82 (9.17)0.80 (0.23)16.66 (21.07)66.94 (19.25)0.20 (0.16)

Table B12 Sample size, mean and variance

Approachnei¯ese2
Stabilized2602.4714.97
Unstabilized6520.80

Table B13 Posterior summaries of parameters for p(Y|A,D,F,G,H,I) with standard deviations (in parentheses)

  Rwy-01Rwy-02Rwy-03
 αy|data3107.63 (53.48)3655.32 (84.21)3553.83 (56.93)
Dβ1y|data10.04 (4.32)110.38 (71.65)91.68 (57.08)
Fβ2y|data5.25 (2.41)3.65 (1.24)4.99 (3.17)
Glβ3y,l|data125.46 (21.04)186.97 (28.57)171.28 (37.80)
Gmβ3y,m|data428.88 (70.88)518.36 (81.12)556.88 (94.58)
Hβ4y|data3.72 (1.31)8.74 (2.33)0.84 (0.56)
Iβ5y|data25.34(12.79)24.11(11.99)23.16(11.06)