A Broad-Based Decision-Making Procedure for Runway Friction Decay Analysis in Maintenance Operations

Abstract

The evaluation of friction is a key factor in monitoring and controlling runway surface characteristics. For this reason, specific airport management and maintenance are required to continuously monitor the performance characteristics needed to guarantee an adequate level of safety and functionality. In this regard, the authors conducted years of experimental surveys at airports including Lamezia Terme International Airport. The surveys aimed to monitor air traffic, features of geometric infrastructure, the typological and physical/mechanical characteristics of pavement layers, and runway maintenance planning. The main objective of this study was to calibrate specific models to examine the evolution of friction decay on runways in relation to traffic loads. The reliability of the models was demonstrated in the light of the significance of the friction measurement patterns by learning algorithms and considering the traffic data by varying the geometric and performance characteristics of the aircraft. The calibrated models can be implemented into pavement management systems to predict runway friction degradation, based on aircraft loads during the lifetime of the surface layers of the pavement. It is thus possible to schedule the maintenance activities necessary to ensure the safety of landing and takeoff maneuvers.

1. Introduction

As a result of the fast growth of the civil aviation industry over the last decades, there have been numerous accidents involving aircraft overruns or deviations from wet or contaminated runways during takeoff and landing. Aircraft operational performance at landing or takeoff is highly dependent on runway surface conditions. For obvious safety reasons, methods and means must be implemented to characterize runway surface conditions and give pilots the relevant information on how well the surface will perform. A plane may find itself off the runway because it has landed at too high a speed or has touched down too late. It is also possible to lose wheel alignment during takeoff and landing. Planes can also lose control as they brake when landing. Runway conditions are a dominant factor in overrun accidents/incidents as the plane comes into land or when taking off.

Runway length has a particular impact on performance during takeoff and landing. Given the shortage of unused space in numerous airports, many are unable to extend their runways, which means that a runway’s friction coefficient and airport clearance are fundamental hardware requirements for airports [1].

Many factors affect friction, including aircraft velocity, wheel load during ground-roll upon landing, the thickness of the water film on the surface, and the pavement surface features. Each plane needs a specific amount of tire–pavement friction for retardation and directional control while taxiing, taking off, and landing. The level of friction required largely depends on the kind of vehicle, its gross weight, the length of the runway, if there is any crosswind, reversed engine thrust usage, and the ability of the pilot [2].

The Netherlands Transport Safety Institute (NLR) carried out research to study the main factors underlying runway excursion accidents in Europe and across the globe [3]. They found that wet or contaminated runways were the leading factor in accidents, making up approximately 40% of European overrun accidents and 75% at the global level. Aircraft braking distance is the prime factor in determining the overall landing distance.

The level of friction available at the aircraft’s tire–pavement interface is particularly important. Substantially reduced pavement skid resistance can be observed in wet weather conditions [4,5].

Furthermore, all the following have been shown to affect wet-pavement skid resistance: tire structure, the characteristics of the pavement surface, the speed of the vehicle, wheel load, tire pressure, and the thickness of the surface water film on the pavement [6,7].

Pasindu et al. [8] developed an analytical simulation model to evaluate the wet-pavement skid resistance available to an aircraft under a given operating condition. The operating condition is defined by aircraft speed, tire structural properties, pavement surface properties, wheel load, tire inflation pressure, and pavement surface water film thickness. The results are used in a numerical example to demonstrate the application of the procedure to evaluate skid resistance variation on a runway during a commercial jet aircraft landing.

Numerous studies have been carried out on tire–pavement friction modeling, adopting approaches that range from semi-empirical methods with experimental curve-fitting to extremely intricate contact models to denote physical occurrences. The frictional features of the interface between the tire and the road are contingent on all the characteristics of the tire (the “rubber texture”, tire geometry, inflation pressure, tread depth, and characteristics of the rubber), operating conditions (slip rate and load), and pavement surface conditions (such as wet, dry, ice, snow, and rubber deposits).

To contribute to tire–pavement friction modeling, Kane and Cerezo [9] proposed an approach based on simplified viscoelastic–rough contact to model the hysteric component of the friction produced by pieces of rubber coming together while transiting on a coarse surface. The rubber’s viscoelastic behavior is modeled using Kelvin–Voigt elements, and the friction force is determined on the basis of the dissymmetry of the contact pressure against asperities on the pavement. A set of friction and texture data was measured on a selection of actual road surfaces in order to confirm the model. Rubber friction has long been an object of study, in particular from the perspective of tire–road interaction. So far, many experimental papers have examined the role of small-scale roughness on rubber friction. However, the possible influences of macro-level texture on rubber friction have often been overlooked. When performing rubber friction experiments on asphalt roads, even if sliding friction experiments are carried out on surfaces with different macro-level textures, it is difficult to attribute the observed friction difference to their macro-texture, because the roughness at other length scales usually differs too.

Kanafin et al. [10] measured friction in relation to two rubber compounds as they slide on rough surfaces randomly produced using 3D-printing technology. The viscoelastic modulus master curves of the rubber compounds and their large strain effective modulus were calculated by means of dynamic mechanical analysis. Numerous prediction models for the skid resistance of asphalt pavements have been developed based on three standard regulation polishing methods: Polished Stone Value, Micro-Deval, and the Wehner/Schulze machine. However, these tests can only be performed under certain conditions, such as constant temperature. This does not wholly represent true road surface conditions because the performance of the asphalt depends very much on temperature on account of the viscosity of the bitumen [11,12].

Xie et al. [13] created prediction skid resistance models and determined the parameters of the models using regression analysis. The authors studied the impact of temperature during polishing on the long-term development of skid resistance. When asphalt surfaces are new, a film of bitumen binder masks the microtexture of the aggregates, which tends to reduce skid resistance in the early life of the pavement. The traffic gradually polishes the microtexture and tends to cause a decrease in skid resistance, which means that pavements will become slippery again with the passage of time.

Kane and Edmondson [14] investigated the possibility of predicting the whole-life-cycle skid resistance of asphalt surfacing. A “weight factor” was added to the dynamic friction model to account for the masking of the aggregate microtexture by a film of bitumen during the early life of an asphalt pavement.

During wintertime, contamination of runway surfaces with snow, ice, or slush causes potential economic and safety threats for the aviation industry. The presence of these materials reduces the available tire–pavement friction necessary for retardation and directional control [15].

The role of airport authorities in guaranteeing that the performance of airport pavements reach the standards required for safe aircraft operations is fundamental [16].

The International Civil Aviation Organization (ICAO) [17] established that the friction value should be obtained by averaging the results of measurements made using the test device. If the friction characteristics differ significantly along major portions of the runway, the friction value should be obtained for each portion. A portion of runway approximately 100 m long may be considered sufficient to establish the friction value. This broad decisional margin in the measurement-processing step is reflected in the identification of the significant value to use in a specific moment of the pavement’s lifetime. Consequently, further knowledge of preliminary data processing methods is required. Airport pavements must ensure regular and safe operations for aircraft. Airport management companies, therefore, have to adopt a method to monitor these surfaces and carry out expensive maintenance and rehabilitation works.

An airport pavement management system (APMS) is an approach to monitoring the condition of pavements and determine the priorities for intervention, as well as planning and allocating resources. An APMS requires multiple tests to evaluate various indices and different skills to conduct appropriate analyses [18].

Airports have dedicated operators and budgets to manage and maintain their airside assets; however, these resources may be lacking in small airports. For managers to prevent the rapid deterioration of the pavement network, it is necessary to implement proactive pavement conservation, applying inexpensive and regular interventions, in contrast with traditional conservation which requires more demanding operations to render pavements useable again [19].

Zou and Madanat [20] examined the way in which an APMS has to consider the scheduling over time of maintenance and rehabilitation actions, and the long-term trade-off between the costs of delayed maintenance and the benefits for air traffic borne by the airport.

De Souza and de Almeida Filho [21] also show the importance of the time factor, suggesting the concept of time delay in maintenance (DTM) as a way of reducing maintenance costs. Their study of a Brazilian airport may be adopted in other contexts after parameterization on the data for failures and runway consumption.

Beyond a certain level, the surface characteristics of a runway can not only affect the costs for upkeep of the pavement, but also that of aircraft landing gear caused by fatigue due to dynamic loads. It is, therefore, essential to find methods able to determine the extent of irregularities (Loprencipe et al. [22]).

For road pavements in general, Pantuso et al. [23] propose modeling a decay index that considers the uncertainties that arise from the acquisition of data and the evaluation of conditions. Decay index curves combine the linear empirical Bayes approach with the measured conditions.

Piryonesi and El-Diraby [24] use causal models to predict pavement condition index (PCI) decay at 2-6 years by using automatic data learning algorithms, two basic decision trees, and their boosted gradient models. With different combinations of 15 attributes concerning the long-term performance of the pavement (LTPP), the gradient boosted models obtained degrees of accuracy greater than 80%.

As for airport infrastructures, Barua et al. [25] use the gradient boosting machine (GBM) learning method in a case study featuring Chicago O’Hare International Airport. They develop two GBM models to study the influence of some attributes on the PCI decay of runway and taxiway pavements: age, type of material, history of maintenance and rehabilitation history, weathering, and air traffic loads. The two models outperformed other methods, both in classical statistics such as linear and nonlinear regressions, and in learning methods such as neural networks and random forests.

Ansarilari and Golroo [26] also use pavement age, maintenance actions, weathering, and traffic load. In their APMS, they use the Markov chain method to model pavement decay. The pavement condition provided by the decay model is then used as an objective function, together with financial resources, to plan maintenance through supervised multi-objective genetic algorithms.

Machine learning algorithms show versatility in many engineering areas that need to extract complex dynamic patterns from time series data. The applications proved to be reliable both regarding issues of predicting the deterioration of individual technological components such as railway track turnouts [27]. They are also reliable in risk-based inspection screening assessment for production and processing units [28] and infrastructure risk assessment in relation to large-scale natural disasters, characterized by data overabundance [29].

Muselli [30] compared two machine learning algorithms, decision trees (DTs) and Hamming clustering (HC), to produce the best reliability expression (RE) estimate.

Hu et al. [31] developed an ensemble of data-driven prognostic algorithms for solid prediction of the useful life remaining, which combine multiple-member algorithms with three weighting-sum schemes: accuracy-based weighting, diversity-based weighting, and optimization-based weighting. Machine learning techniques are also used in non-engineering systems, such as companies assessing corporate governance.

Hernandez-Perdomo et al. [32] sampled 1109 companies listed on the United States of America stock exchange from 2002 up to 2014, setting up the components of the companies and the company characteristics as input and the economic performance as output. A structure function was, therefore, defined to model the performance of the activities.

The authors have been conducting years of experimentation surveys at airports [33,34,35]. The surveys aimed to monitor air traffic, features of geometric infrastructure, the typological and physical/mechanical characteristics of pavement layers, and runway maintenance planning.

The aim of this study was to carry out a theoretical–experimental analysis of the evolution of friction decay on runways, focusing on developing procedures for the preliminary analysis of the raw data collected and on the construction of models to reliably represent the evolution of the surface friction characteristics on the basis of traffic loads. The reliability of the friction models was demonstrated in the light of the significance of the friction measurement patterns by learning algorithms: a first part of the analysis looked at the raw friction data for which different pattern recognition methods were selected and tested to locate the consumption zones with their spread and intensity. The reliability of the traffic data was tested by varying the geometric and performance characteristics of the aircraft: this part of the analysis looked at the aircraft traffic mix in terms of configurations, flight plans, and engine performance to calculate actual transferred loads to the pavement. The calibrated models may be implemented into pavement management systems to predict runway friction degradation based on aircraft loads during the lifetime of the surface layers of the pavement. It is thus possible to schedule the maintenance activities necessary to ensure the safety of landing and takeoff maneuvers.

2. Methods and Analysis

2.1. Data Processing Methods

Airport pavement management systems (APMSs) can only benefit from data if the findings are studied in such a way that they provide a reliable image of how decay develops. The principal tools to map decay are the grip number (GN) and time decay curves, these need a fixed value for every relevant date. From the analysis of the available ICAO standards [17,36,37,38], it can be assumed there are no specific indications on the exact section lengths of the runway for the assessment of friction or on the methods of choosing the measurement to which the level of degradation can be referred.

To evaluate the best method to establish this measurement value, the authors selected four data mining methods, namely computational processes for pattern recognition in large data sets. Two supervised learning algorithms were applied: classification and regression trees (CARTs) [39] and general chi-squared automatic interaction detector (GCHAID) [40]; and two unsupervised learning algorithms: k-means [41,42,43,44] and Chiu’s subtractive clustering [45,46]. The results of the numerical processing were compared with what could be obtained by applying the methods commonly used in technical practice for the same purpose as the one-third segment (ENAC APT-10A) [47] and the Minimum 100 Meters Rolling Averages (CAA M100MRA) [48] methods.

Then, air traffic data were processed in order to obtain a friction decay model based on actual weights (instead of nameplate weights): landing and takeoff weights were estimated on flights over 12 years, and the historical series of these cumulated weights were analyzed.

Based on actual landing weights and on geometric configuration and performance characteristic of the aircraft spectrum, a model for cumulated load for each pavement alignment was calibrated. Then, an empiric-mechanistic friction-load model was optimized and validated, using friction data sample excluded from the calibration phase, through residuals analysis.

The last part of the process was an application of the models for airport pavement maintenance purposes: assessment of the maximum cumulated load alignment, calculation of the cumulative load on that alignment, and scheduling the maintenance operations.

The methodological process is shown in Figure 1.

2.2. Case Study General Overview

The case study was Lamezia Terme International Airport, a medium-sized airport with a passenger count of 2,756,211 recorded for 19,102 flights in 2018 [49]. The runway, named 10/28, based on the magnetic orientation criterion, is classed as a 4D on the ICAO scale, with a field length of over 1800 m, accommodating wingspans ranging from 36 to 52 m and a main landing gear span between 9 and 14 m. Air traffic may follow both instrument flight rules (IFR) and visual flight rules (VFR).

The runway has a total paved area of 60 m × 2414 m = 144,840 m²; the part bounded by the runway thresholds, which are the markings across the runway that denote the beginning and end of the designated space for landing and takeoff under non-emergency conditions, is 45 m wide by 2200 m in length.

The pavement classification number, which is a number expressing the bearing strength of a pavement for unrestricted operations, is PCN 95/F/B/W/T.

The pavement surface layers were rebuilt in 2004 during general refurbishment work on most of the airside infrastructure. After the pavement was treated, a 12 cm scarification was made, with a 12 cm refill using a bituminous conglomerate as binder and a 5 cm weak layer containing modified bitumen on a 25-meter band straddling the runway’s axis line. A 5 cm scarification was made on the adjacent 10-meter strips, with a 5 cm refill using a bituminous conglomerate as binder and a 5 cm weak layer with modified bitumen.

The measurements of the grip number (GN) were made along the alignments at ±3 m and ±6 m from the centerline, which are the most used alignments reflecting the main landing gears of the aircraft traffic mix, and for the whole length of the runway. Data were collected using a Grip Tester Trailer device with a 10 m spacing output [17]. The final database consisted of 92 alignments of raw data, which refer to 23 surveys made from 2004 to 2015 (see Supplementary Materials, S1).

In July 2007 and November 2011, derubberizing maintenance operations were carried out using a high-pressure water TrackJet machine.

Constant GN values were assigned for 10 and 28 runway thresholds.

In order to eliminate the effects of the transient speed of the Grip Tester at start (before one threshold) and at stop (after the other threshold), those GNs were set at the first value on the respective threshold. For this reason, only the values of the runway area bounded by the horizontal marking of the two thresholds were processed.

2.3. Friction Data Analysis

The main objective in this phase was to segment the runway longitudinally into parts reflecting the diverse zones of decay on each survey date in order to identify the one with the highest level of degradation.

Ninety-two data series based on 23 surveys for the alignments at ±3 m and ±6 m were collected, of which 76 were used in the calibration phase and the remaining 16 in the validation phase.

In the processing steps, the friction values for each of the two symmetrical alignments, +3 and −3, +6 and −6, were averaged out and will hereafter be referred to as ±3 and ±6.

Table 1. Overview of the GN minimum values.

Day	CART		GCHAID		K-means		Subtractive Clustering		APT-10A		M100MRA
	±3 m	±6 m	±3 m	±6 m	±3 m	±6 m	±3 m	±6 m	±3 m	±6 m	±3 m	±6 m
03/08/2004	-	0.90	-	0.91	-	0.90	-	0.90	-	0.17	-	0.89
08/31/2004	-	0.88	-	0.94	-	0.87	-	0.93	-	0.67	-	0.87
06/16/2005	-	0.86	-	0.89	-	0.81	-	0.90	-	1.46	-	0.84
12/21/2005	-	0.85	-	0.91	-	0.85	-	0.90	-	2.00	-	0.81
05/31/2006	-	0.81	-	0.83	-	0.79	-	0.85	-	2.42	-	0.78
07/03/2007	-	0.72	-	0.74	-	0.71	-	0.68	-	3.50	-	0.69
10/24/2008	0.82	0.85	0.93	0.99	0.81	0.86	0.84	0.87	0.94	4.83	0.78	0.82
07/22/2009	0.75	0.81	0.85	0.94	0.75	0.80	0.69	0.82	0.86	5.58	0.68	0.77
10/14/2009	0.69	0.77	0.80	0.90	0.69	0.76	0.66	0.77	0.81	5.83	0.63	0.72
10/28/2010	0.57	0.72	0.74	0.87	0.56	0.72	0.63	0.74	0.75	6.83	0.57	0.69
04/29/2011	0.50	0.69	0.62	0.84	0.47	0.68	0.51	0.71	0.64	7.33	0.44	0.66
05/27/2011	0.43	0.64	0.58	0.77	0.45	0.61	0.46	0.65	0.58	7.41	0.38	0.60
10/11/2011	0.38	0.60	0.47	0.74	0.38	0.60	0.35	0.60	0.47	7.83	0.36	0.57
11/15/2011	0.33	0.56	0.40	0.62	0.34	0.57	0.29	0.54	0.41	7.91	0.29	0.52
02/15/2012	0.73	0.81	0.71	0.83	0.73	0.82	0.77	0.85	0.76	8.16	0.71	0.80
05/17/2012	0.72	0.80	0.73	0.84	0.75	0.84	0.74	0.87	0.74	8.41	0.70	0.79
10/24/2012	0.69	0.79	0.70	0.82	0.70	0.78	0.73	0.83	0.73	8.83	0.68	0.78
12/21/2012	0.68	0.78	0.72	0.80	0.69	0.77	0.70	0.83	0.72	9.00	0.67	0.77
04/18/2013	0.67	0.77	0.72	0.80	0.71	0.77	0.72	0.79	0.71	9.33	0.66	0.76
05/16/2013	0.64	0.74	0.68	0.79	0.61	0.74	0.71	0.72	0.69	9.41	0.61	0.71
12/16/2013	0.59	0.68	0.71	0.78	0.57	0.69	0.60	0.69	0.68	10.00	0.56	0.68
05/28/2014	0.51	0.63	0.61	0.72	0.51	0.64	0.64	0.64	0.62	10.42	0.48	0.59
02/11/2015	-	-	-	-	-	-	-	-	-	-	-	-

shows an overview of the GN minimum values obtained by varying the algorithms used for each survey and alignment (±3 m; ± 6 m). The empty fields, marked “-”, refer to the 7 data series used in the validation phase, namely the measurements collected for the ±3 m alignments from 2004 to 2006 and the ±3 m/±6 m alignments during the last survey in 2015.

The GN values, obtained by applying the 6 methods listed in the methods section, were used to construct time–GN decay curves for alignment ±3 m (see Figure 2a) and alignment ±6 m (see Figure 2b). The second-degree polynomial function proved to be the most suitable for all the diagrams, generating higher R² determination coefficients than other linear or nonlinear functions.

A general shape shown by all curves is due to maintenance operations carried out in July 2007 and November 2011: after the derubberizing maintenance operations, there was a sharp increase in GN, dividing each curve into separate branches (2 branches for the ±3 m alignment and 3 branches for the ±6 m alignment).

By comparing the time–GN decay curves obtained using the different methods, it was possible to make the following observations:

• CART decay curves show regular and comparable slope patterns in both ± 3 m and ± 6 m branches. The GN values are always conservative (i.e., low) and progressive.
• GCHAID decay curves show fewer regular patterns. The 2008–2011 ±6 m branch starts higher than the two adjacent ones (2004–2007 and 2012–2014). The 2012–2014 ±3 m and ±6 m branches have significantly lower slopes than previous branches. The GN values are less conservative and, in some cases, stationary and recoiling.
• K-means decay curves have similar slope patterns to CART decay curves and conservative GN values, and in some cases, they are stationary and recoiling.
• The subtractive clustering decay curves appear to be less regular than the CART and K-means decay curves. The 2012–2014 ±3 m branch has a nearly constant slope (a linear trend), and the 2012-2014 ±6 m branch appears to be more upwards moving than the previous ones. The GN values are in some cases stationary or recoiling.
• APT-10A decay curves show irregularities like the GCHAID and K-means decay curves. The GN values are always progressive, but they are the least conservative (i.e., the highest), being the result of averages on very large parts of the runway.
• M100MRA decay curves appear regular, though less than CART and K-means decay curves, especially at alignment ±6 m. They have progressive values.

Four criteria were established in order to carry out a quantitative comparison of the previous observations, where scores or marks were assigned on a scale from 0 to 10, and weights were assigned on the basis of interviews with researchers and professional technicians operating in the field of airport infrastructures. They evaluated the different criteria in terms of their relative importance.

In particular:

• The conservative nature of terminal values: a decay curve is as reliable as the GN values identified by its lower algorithm and is, in this way, closer to the scheduled minimum maintenance levels set by the standards for the specific means of measurement. The maximum score is 10 and is attributed to the decay curve with the lower mean terminal value and by correlating the others according to their mean terminal values. Weight value P is 3.
• The shape of the decay curve: a decay curve is reliable if its branches have a regular layout and slope. The maximum score, essentially qualitative, is 10, assigned to the CART decay curves and by deducting one point for each identified irregularity. Weight value P is 1.
• The progression of values: a decay curve is reliable if its values decrease over time, with no stationary or recoiling values as time goes by. The maximum score is 10 for the decay curves with all its key values in progression, deducting one point for each stationary value and/or instance of recoiling. Weight value P is to 2.
• Standard deviation distribution: a decay curve is reliable if the single point values from which it is derived (for non-linear regression) are accurate, that is with a low standard deviation. The maximum score is 10 for the curve of each alignment with the lowest standard deviation distribution and taking the others for each alignment according to the respective Standard deviation distribution, having chosen the average of distribution as a parameter. Weight value P is 2.

Table 2. Summary of the multicriteria analysis.

Criteria		The Cautiousness of Terminal Values	Shape of Decay Curve	Progression of Values	Standard Deviation Distribution	Standardized Score
CART	±3 m	9.39	10	10	9.20	12.44
	±6 m	9.57	10	10	8.92	12.09
GCHAID	±3 m	7.62	9	6	4.58	8.16
	±6 m	8.66	7	6	5.02	8.29
K-means	±3 m	9.28	10	8	10.00	12.09
	±6 m	9.52	10	7	10.00	11.55
Subtractive Clustering	±3 m	8.28	9	8	4.67	9.13
	±6 m	9.68	7	6	4.74	8.52
APT-10A	±3 m	7.48	9	10	4.35	9.33
	±6 m	8.45	8	10	4.54	9.48
M100MRA	±3 m	10.00	10	10	5.78	11.11
	±6 m	10.00	9	10	6.16	10.88

shows the results obtained for alignment ±3 m and ±6 m using multicriteria analysis.

For both alignments ±3 m and ±6 m the CART method was more reliable than the other five and was therefore selected for further analysis.

The main output at the end of the CART implementation is the graph of the optimal tree structure, with the terminal nodes labeled with the average of the GN values of the elements included in the node. We set up GN as the dependent variable, the metric progressive of the runway from threshold 10 as the predictor variable, taking a minimum of 22 cases as a stopping rule, together with a 10-fold cross-validation, and the 1-SE pruning rule [39]. Figure 3 shows the optimal tree relating to the series of the 4/29/2011 survey on ±6 m alignments.

The variance of the elements of the node and the number of elements are also within each node. Above each node is the split value of the progressive metric, i.e., the decision rule that partitions the elements of the node based on the value of the distance from threshold 10. This output information is meaningful with respect to the decay phenomenon: it locates the various friction consumption areas along the runway and indicates their extension.

Figure 4 shows the friction consumption zones for the series of the 4/29/2011, 5/27/2011, and 10/11/2011 surveys on ±6 m alignments from the CART algorithm: those with maximum consumption (i.e., with the lowest GN) are labeled with their GN, spread, and localization.

2.4. Traffic Data Analysis

2.4.1. Estimating Actual Weights

Traffic data for the period 2004–2015, provided by airport management company S.A.CAL. S.p.A., were processed in order to obtain a reliable estimate for the actual weights of the aircraft for each landing and takeoff.

These data are not directly available, unlike nameplate data (i.e., maximum landing weight— MLW): nameplate weights are fixed limit weights used in airport design and operation due to their immediate availability. In order to obtain a mechanistic load model, it was decided to introduce an estimate of the actual weights.

The most flexible parameter for calculating these weights is the amount of fuel in the tanks. This parameter was estimated on the basis of flight plans, in terms of fuel consumption, and according to the standard recommendations on the amount of fuel to be insured, taking two observed conditions into account [50]: the average traffic level, and, therefore, the engagement of the runways, is 3 flights per hour, 16 hours a day, which implies almost no “holding” at the specified 1500 ft quota. This is reflected in the lack of consumption of the final reserve fuel and additional fuel of the landing aircraft.

The geographical location had favorable meteorological conditions prevailing during the year, which implied a very low percentage of redirects (which may occur in the nearby airports of Reggio Calabria and Naples). This is reflected in flight plans without the need for an alternative aerodrome.

The landing weights (LWs) were calculated as:

$$\mathrm{LW}=\mathrm{Operational}\mathrm{Empty}\mathrm{Weight}\left(\mathrm{OEW}\right)+\mathrm{Pay}\mathrm{Load}\left(\mathrm{PL}\right)+\mathrm{Landing}\mathrm{Fuel}\left(\mathrm{LF}\right)$$

(1)

where OEW is the weight of the aircraft under operating conditions ready to fly with fluids and crew included, excluding fuel and payload; PL is passengers with luggage estimated for aircraft with 30 or more seats, since Lamezia Terme Airport is passenger-only traffic. The passenger mass value is equal to 76 kg, and for aircraft with 20 or more seats the standard flight baggage for European regions equals 13 kg; LF is reserve fuel, sum of the contingency fuel + final reserve fuel + additional fuel, in practical terms equal to 5% trip fuel + 45 min/30 min fuel + 15 min fuel. Trip fuel is fuel consumption (kg/h) × flight time between the airport of origin and Lamezia Terme Airport (h) at cruising speed. The 45 min/30 min fuel is the fuel for an additional 45 min/30 min flight time for piston/turbine-powered aircraft.

$$\mathrm{TOW}=\mathrm{Operational}\mathrm{Empty}\mathrm{Weight}\left(\mathrm{OEW}\right)+\mathrm{Pay}\mathrm{Load}\left(\mathrm{PL}\right)+\mathrm{TakeOff}\mathrm{Fuel}\left(\mathrm{TOF}\right)$$

(2)

In the case of TOF, it is not possible to assume the same hypothesis as for LF since destination airports may present one (or even two) alternative airports, depending on their traffic characteristics and weather conditions. The simplifying hypothesis TOF = trip fuel + reserve fuel = additional fuel = trip fuel + 5% trip fuel + 45 min/30 min fuel + 15 min fuel will be used, but as will be seen later, this simplification will not affect the mechanistic load model.

The estimates of real weights for 72,572 landings and takeoffs were calculated for the period 2004–2015.

S2 in Supplementary Materials includes the nameplate data, configurations, and performance data for the spectrum of the 72 aircraft at Lamezia Terme Airport from 2004 to 2015, while S3 in Supplementary Materials, an excerpt of the traffic data calculated using the estimates of the above weights in [kg], refers to the period between the surveys of 7/22/2009 and 10/14/2009.

2.4.2. Historical Series of Cumulated Loads and Identification of the Regressive Variable

Cumulative landing weight (CLW) and cumulative takeoff weight (CTOW) were then calculated for the period 2004–2015. The Pearson coefficient calculated using Equation (3) (where: x_i = i—cumulative landing weight; y_i = i—cumulative takeoff weight) confirmed a very significant correlation between CLW and CTOW (ρ_xy = 0.999).

$$-1\le {\mathsf{\rho }}_{\mathrm{x}\mathrm{y}}=\frac{{\mathsf{\sigma }}_{\mathrm{x}\mathrm{y}}}{{\mathsf{\sigma }}_{\mathrm{x}}{\mathsf{\sigma }}_{\mathrm{y}}}=\frac{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}({\mathrm{x}}_{\mathrm{i}}-{\mathsf{\mu }}_{\mathrm{x}})({\mathrm{y}}_{\mathrm{i}}-{\mathsf{\mu }}_{\mathrm{y}})}}{\sqrt{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{x}}_{\mathrm{i}}-{\mathsf{\mu }}_{\mathrm{x}})}^2}}\sqrt{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{y}}_{\mathrm{i}}-{\mathsf{\mu }}_{\mathrm{y}})}^2}}}\le 1$$

(3)

The independent variable between CLW and CTOW in the model was decided based on mechanics-related considerations. Given the two rolling conditions of the wheels during landing (braked wheel) and takeoff (towed wheel), in the case of a braked wheel, friction force A expressed at tire–pavement contact equals to:

$$\mathrm{A}=\left({\mathrm{M}}_{\mathrm{a}}+{\mathrm{M}}_{\mathrm{f}}\right)/\mathrm{r}$$

(4)

where, M_a is the frictional moment at the wheel axle pin, M_f is the braking moment, and r is the radius of the wheel; while in the case of a towed wheel it is always equal to:

$$\mathrm{A}={\mathrm{M}}_{\mathrm{a}}/\mathrm{r}$$

(5)

As M_f is much greater than M_a, the friction forces that contribute significantly to the degradation of the surface characteristics of the pavement are explained during the landing phase. Based on the above considerations, CLW was considered an independent variable of the mechanistic load model to be calibrated. Twelve different aircraft models made up 93.3% of the flights in the study period as follows: Boeing B737-700 (B737), McDonnel Douglas MD82 (M82), McDonnel Douglas MD80 (M80), Airbus 320 (A320), Boeing B737-800 (B738), Airbus 319 (A319), Fokker 100 (F100), Airbus 321 (A321), Boeing B737-400 (B734), Bombardier CRJ900 (CRJ), Dornier 328 (D328), Tupolev 154 B2 (TU154). For each of these, a comparison was carried out between LW and MLW, returning LW values lower than MLW in every case, with very pronounced seasonal trends.

Table 3. Summary of the main characteristics of the 12 aircraft types.

	No. Flights	LWm	s	C.V.	MLW	MLW-LWm	(MLW-LWm)/LWm
B737	16,652	48,127	3970	0.082	58,060	9933	0.206
M82	12,637	47,386	3692	0.078	58,970	11,584	0.244
M80	10,811	45,405	4891	0.108	58,060	12,655	0.279
A320	9655	53,262	3685	0.069	64,500	11,238	0.211
B738	5891	57,334	3416	0.060	65,310	7976	0.139
A319	4460	51,930	2729	0.053	61,000	9070	0.175
F100	3023	33,293	2540	0.076	38,780	5487	0.165
A321	2343	62,341	4366	0.070	75,500	13,159	0.211
B734	1257	45,902	3889	0.085	54,880	8978	0.196
CRJ	651	27,733	1920	0.069	33,339	5606	0.202
D328	272	11,698	934	0.080	14,390	2692	0.230
TU154	63	64,790	6889	0.106	80,000	15,210	0.235

Abbreviations: C.V.—Coefficient of Variation, LW—landing weight, MLW—maximum landing weight, B737—Boeing B737-700, M82—McDonnel Douglas MD82, M80—McDonnel Douglas MD80, A320—Airbus 320, B738—Boeing B737-800, A319—Airbus 319, F100—Fokker 100, A321—Airbus 321, B734—Boeing B737-400, CRJ—Bombardier CRJ900, D328—Dornier 328, TU154—Tupolev 154 B2.

shows a summary of the main characteristics of the 12 aircraft types. In particular, it was observed that the average difference between the MLW and the mean value of the LW (LW_m) is 20.8%, and the absolute differences between MLW and LW_m increase with landing capacity depending on the MLW. This confirms the reliability of using LW in the modeling phase.

2.4.3. The Mechanistic Load Model with Actual Weights and Calibration

The mechanistic loads model was derived from the LW using Equation (6). The cumulated load of j-aircraft on each measurement alignment y takes into account the geometric configuration of the landing gears, the cross dispersion of the landing trajectories, the distance between the runway centerline, and the measurement alignments.

For the configuration of the main landing gear (MLG) in 72,033 of the 72,572 landings considered, i.e., with a two-strength leg (right and left form the aircraft axle) and wheels on two alignments for each strength leg (outer wheels and inner wheels):

$${\mathrm{CLW}}_{\mathrm{y}}={\displaystyle \sum _{\mathrm{j}}{\mathrm{LW}}_{\mathrm{jy}}=}{\displaystyle \sum _{\mathrm{j}}{\mathrm{LW}}_{\mathrm{j}}·0.95·0.25·{\displaystyle \sum _{\mathrm{k}=\mathrm{outer}}^{\mathrm{inner}}{\displaystyle \underset{\mathrm{y}-\frac{\mathrm{d}}{2}}{\overset{\mathrm{y}+\frac{\mathrm{d}}{2}}{\int }}\mathrm{pdf}\left[{\mathsf{\mu }}_{\mathrm{j}}=\pm {\mathrm{Wa}}_{\mathrm{kj}}{\mathrm{Aa}}_{\mathrm{j}};\mathsf{\sigma }\right]}}}$$

(6)

where, LW_jy is the landing weight of the j-aircraft on the measurement alignment y; LW_j is the landing weight of the j-aircraft; 0.95 is the LW_j distribution coefficient on the MLG; 0.25 is the distribution coefficient on the MLG’s single wheel alignment (right outer, right inner, left inner, or left outer); y is the distance from the centerline of the measurement alignment (3 m or 6 m); d is the thickness of the measuring wheel of the Grip Tester Trailer (equal to 4 inch = 0.1016 m); pdf is the Gaussian probability density function; μ is the Gaussian parameter; Wa_kjAa_j is the difference between the aircraft vertical axle and MLG wheel vertical axle (placing the origin of the crossing distances from the centerline, as we averaged out the friction values for the symmetrical alignments, +Wa_outerjAa_j and +Wa_innerjAa_j are for the outer and inner wheels of the right strength leg, -Wa_innerjAa_j and -Wa_outerjAa_j are for the inner and outer wheels of the left strength); σ is the Gaussian standard deviation of the landing trajectories across the entire traffic spectrum (for the lateral wander effect).

Equation (7) can also be used:

$$\mathrm{CL}{\mathrm{W}}_{\mathrm{y}}={\displaystyle \sum _{\mathrm{j}}{\mathrm{LW}}_{\mathrm{jy}}=}{\displaystyle \sum _{\mathrm{j}}\mathrm{L}{\mathrm{W}}_{\mathrm{j}}·0.95·0.5·{\displaystyle \underset{\mathrm{y}-\frac{\mathrm{d}}{2}}{\overset{\mathrm{y}+\frac{\mathrm{d}}{2}}{\int }}\mathrm{p}\mathrm{d}\mathrm{f}\left[{\mathsf{\mu }}_{\mathrm{j}}=\pm \mathrm{S}{\mathrm{a}}_{\mathrm{j}}\mathrm{A}{\mathrm{a}}_{\mathrm{j}};\mathsf{\sigma }\right]}}$$

(7)

where, 0.5 is the distribution coefficient on the MLG’s single-strength leg (right or left); Sa_jAa_j is the difference between the aircraft vertical axle and the MLG axle leg strength (+Sa_jAa_j for the right strength leg, -Sa_jAa_j for the left strength leg). The calculations showed imply differences between LW_jy by Equation (6) and LW_jy by Equation (7) of −0.003177 = −0.32% on the alignment +3 m and of 0.004946 = 0.49% on the alignment +6 m.

Equation (7) identifies the load for MLG axle leg strength as responsible for friction consumption. Once the weights and geometric configuration of the aircraft are established, Equation (7) returns CLW_y for the alignments at ±3 m and ±6 m, varying on the basis of σ. The σ parameter was calibrated by maximizing the index of determination of the regression line of the initial points CLW_y(σ)–GN of each decay curve branch ±3 m and ±6 m (see Figure 5). The starting points of the branches represent the post-derubberizing initial states. This is due to the physical behavior of the phenomenon, confirmed by the data: post-derubberizing returns the surface elements to a "renewed" state whose level, that is the extent of the recovery in GN, depends on the history of the loads and not on the alignment.

The regression for the initial points of the branches is linear, being the most suitable for the 2 initial points of the time–GN ±3 m decay curve and the 3 initial points of the time–GN ±6 m decay curve. CLW_y–GN decay curves can be plotted on the same diagram without distinction of alignment. The optimal value for σ is 3.411 m.

2.4.4. The Empirical-Mechanistic Friction Decay Model

A non-linear regression analysis makes it possible to obtain the mathematical function of the 5 branches of the 2 experimental CLW_y–GN decay curves referring to ±3 m and ±6 m. The second-degree polynomial function proved to be the most suitable for all the branches, generating determination coefficients R² higher than other linear or nonlinear functions. The subsequent optimization procedure obtained the variability laws of the parameters a₁, a₂, and a₃ of the second-degree polynomial functions on the CLW_y–GN plane.

Here is the empirical-mechanistic friction model, Equation (8):

$$\mathrm{GN}={\mathrm{GN}}_{\mathrm{new}}-\left(\Delta {\mathrm{GN}}^{\prime }+\Delta {\mathrm{GN}}^{″}\right)={\mathrm{GN}}_{\mathrm{restart}}-\Delta {\mathrm{GN}}^{″}$$

(8)

where, GN_new is the GN value for the new pavement; GN_restart = −0.0121CLW_y + 0.905 is the GN value after derubberization; ∆GN’ is the loss of non-recoverable friction calculated as the difference between GN_new and GN_restart; ∆GN’’ is the loss of recoverable friction calculated as shown in Equation (9):

$$\Delta \mathrm{G}{\mathrm{N}}^{″}={\mathrm{a}}_1{\mathrm{CLW}}_{\mathrm{y}}^2+{\mathrm{a}}_2{\mathrm{CLW}}_{\mathrm{y}}+{\mathrm{a}}_3$$

(9)

where, a₁ = −0.41905GN_restart² + 0.6179GN_restart − 0.2299; a₂ = −82.445GN_restart³ + 196.49GN_restart² − 155.64GN_restart + 41.03; a₃ = −4034GN_restart⁴ + 13281GN_restart³ − 16351GN_restart² + 8923.8GN_restart − 1821.35.

The model shows the loss of friction not recoverable through derubberization, attributable to textural consumption, as well as temporary loss of friction recoverable through derubberization, attributable to the texture covered by gummy deposits. The functional form of ∆GN’’ indicates that it depends not only on the accumulation of load but also on the starting points of the curves, namely the initial post-derubberizing state. Different GN values can be obtained for the same CLW_y value on the abscissa: this depends on which point of the previous curve the surface element was when the maintenance operation was performed, and then the regained position on the GN_restart line. In other words, it depends on the history of the accumulated loads and maintenance operations performed.

A graduated abacus on the plane CLW_y–GN was obtained from the empirical-mechanistic model shown in Figure 6. It can account for the phenomenon visually and be used for graphic calculation. The legend to the right shows the restarting points (CLW–GN) of each plotted branch placed on the restarting line.

2.5. Validation Phase

Validation of the empirical-mechanistic model was carried out by comparing the results predicted by Equation (8) with 8 experimental points: 6 deriving from the 2004–2007 surveys on the ±3 m alignment and 2 from the 11/02/2015 survey on both the ±3 m and ±6 m alignments. CLW on ±3 m and ±6 m alignments were calculated through Equation (7), see

Table 4. Predicted GN versus observed GN.

Day	CLW [Kt]		GN Predicted		GN Observed
	±3 m	±6 m	±3 m	±6 m	±3 m	±6 m
03/08/2004	0.210		0.913		0.900
08/31/2004	0.831		0.893		0.870
06/16/2005	1.972		0.827		0.800
12/21/2005	2.905		0.744		0.760
05/31/2006	3.607		0.665		0.630
07/03/2007	5.456		0.388		0.410
02/11/2015	20.789	12.269	0.459	0.599	0.470	0.570

.

The validation procedure returned the following results:

• A performance diagram where the observed GN lies on the x-axis and the predicted GN is on the y-axis; the layout of the points is fairly well aligned to the bisector;
• Assessment of the mean residual GN_predicted–GN_observed at 0.0098 GN and standard deviation at 0.0226;
• Analysis of the standard deviation (σ) and mean (μ) of the residualizing of the “3σ” method; all residuals fall within the range [μ − 2σ; μ + 2σ], where μ is the mean value and σ is the standard deviation of the residuals;
• Assessment of the accuracy indicators MAD (mean absolute deviation), MAPE (mean absolute percentage error), and MSE (mean squared error) at 0.0219, 3.4816, and 0.0005 respectively, confirming the statistical significance of the prediction model.

3. Results

A preliminary statistical analysis of the parameter SajAaj was carried out in order to identify the maximum cumulated load’s alignment. The result showed peaks of absolute frequencies corresponding to MLG with SajAaj equal to 1.78 m (absolute frequency n = 23,594), 2.30 m (absolute frequency n = 22,787) and 3.19 m (absolute frequency n = 16,188). The mean, mode, and median values of the entire distribution were 2.37 m, 1.78 m, and 2.30 m respectively.

CLW ± 1.78, CLW ± 2.30, and CLW ± 2.37 were calculated using Equation (7). The results returned the maximum cumulated load on the alignment at ±2.37 m from the centerline, i.e., at a distance from the centerline equal to the mean of the SajAaj of the MLG of the aircraft.

The historical series of cumulative loads transversal to the runway calculated for the following alignments 0 m, 1.50 m, 2.37 m, 3.00 m, 4.50 m, 6.00 m, 7.50 m, 9.00 m, 12 m, 15 m is shown in Figure 7.

It is possible to schedule a pavement maintenance operation on the basis of the maximum cumulation alignment ±2.37 m. Friction levels represented by the decay curve must be compared with the following friction levels provided in the ENAC (i.e., Italian Civil Aviation Authority) Guide [47]: Maintenance Planning Level (MPL) at 0.53 and Minimum Friction Level (MFL) at 0.43. By always scheduling derubberizing maintenance operations when MFL is reached, it is possible to establish the pair of values CLW_±2.37–GN to obtain the GN_restart value, as shown in Figure 8.

From the historical series of cumulative loads CLW_±2.37 shown in Figure 7, plotted against time in Figure 9, it is possible to infer what should have been the correct derubberizing dates considering the ±2.37 m alignment rather than ±3 m:

• 1^st derubberizing: CLW_±2.37 = 5.2 kt → 3.35 years, i.e., April 2007 (instead of July 2007);
• 2^nd derubberizing: CLW_±2.37 = 11 kt → 6.5 years, i.e., June 2010 (instead of November 2011);
• 3^rd derubberizing: CLW_{± 2.37} = 18.4 kt → 9.95 years, i.e., December 2013 (not yet performed by February 2015).

In the same way, using the equation shown in Figure 9, it is possible to schedule further maintenance operations.

4. Discussion and Conclusions

The main objective of this study was to calibrate specific models to examine the friction decay evolution on runways in relation to traffic loads, focusing on developing procedures for the preliminary analysis of the raw data collected and on the construction of models to reliably represent the evolution of the surface friction characteristics.

From the study of the current regulations, it proved necessary to deepen knowledge of the preliminary procedures for processing the relevant raw data relating to the surface characteristics of runways.

Measurements obtained from high-performance instruments like continuous friction measurement equipment supplied a sizeable database whose data can be mined to provide an accurate description of the physical phenomenon in question.

The method proposed by the authors makes it possible to identify the position and extension of the maximum friction consumption area and to recognize its “label” value as the minimum reference friction value for the entire runway based on requirements and criteria such as the possibility of identifying spatially continuous consumption zones and time-varying localization and extension, conservative label values, the reliability of decay curve layouts, the progression of values, and precision. This preliminary step of analyzing the relevant raw data proved necessary to correctly represent the evolution of the phenomenon at different times in the useful life of the pavement. The analysis showed that data mining using a CART learning algorithm returned particularly consistent and accurate findings. The algorithm made the data self-standing, elastic enough to evaluate the various zones of consumption from the perspective of both spread and intensity.

By analyzing actual loads, it was possible to calculate an empirical-mechanistic model, producing the σ parameter (transversal dispersion of the trajectories in touchdown on the entire traffic spectrum), which is also innovative with regard to the actual value of this parameter.

The proposed methodology was based on solid mathematical models, although newly applied in this field, on standardized and widely used measuring instruments in the airport sector and on the geometric configuration and performance characteristics of aircraft.

Regardless of ICAO class, the proposed methodology can be applied to other medium-sized airports with passenger-only traffic and flexible pavement design. The aircraft spectrum, explored more in-depth in Supplementary Materials, S2, cover a range of 3 to 440 seats; 1035 kg to 330,000 kg for maximum landing weight; 2 to 4 for rear landing gear; 2.30 m to 11.50 m for track; 2;4 to 4;20 for wheels (nose; main); and single to tri-twin-tandem for type of wheels.

In conclusion, the engineering contribution of this model to APMSs consists of its use in the planning of ordinary and extraordinary maintenance activities, necessary to ensure safety. Once the maximum cumulation alignment has been identified on the basis of traffic data, the friction decay phenomenon and the correct dates on which maintenance operations can be better forecasted.

Another reliable result is the validation of the cumulative character of the friction decay phenomenon, making it possible to distinguish between a temporary and recoverable loss of friction due to the deposit from rubber tire contaminants and a non-recoverable loss of friction caused by the degradation of the surface texture.

Given the proposed methodology, downstream of this study, future development emerges in validating the models for other runways. How to derive and validate other estimates of actual landing weights is also a possible object of new study, for example by properly evaluating fixed plate data such as maximum landing weights, or by other projections based on similar traffic data from other airports.