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

: The evaluation of friction is a key factor in monitoring and controlling runway surface characteristics. For this reason, speciﬁc 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 tra ﬃ c, 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 speciﬁc models to examine the evolution of friction decay on runways in relation to tra ﬃ c loads. The reliability of the models was demonstrated in the light of the signiﬁcance of the friction measurement patterns by learning algorithms and considering the tra ﬃ c 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 takeo ﬀ maneuvers.


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.

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. 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 and then 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.

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

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 2 ; 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.

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 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 2 determination coefficients than other linear or nonlinear functions.
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 2 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).
(a) 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 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 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. Sustainability 2020, 12, x FOR PEER REVIEW 10 of 21 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.

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. 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.  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.

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.

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: 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 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: 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: 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.  Table 3 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.

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): 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 kj Aa 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 outerj Aa j and +Wa innerj Aa j are for the outer and inner wheels of the right strength leg, -Wa innerj Aa j and -Wa outerj Aa 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: where, 0.5 is the distribution coefficient on the MLG's single-strength leg (right or left); Sa j Aa j is the difference between the aircraft vertical axle and the MLG axle leg strength (+Sa j Aa j for the right strength leg, -Sa j Aa 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. Sustainability 2020, 12, x FOR PEER REVIEW 14 of 21 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. CLWy-GN decay curves can be plotted on the same diagram without distinction of alignment. The optimal value for σ is 3.411 m.

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 CLWy-GN decay curves referring to ±3 m and ±6 m. The seconddegree polynomial function proved to be the most suitable for all the branches, generating determination coefficients R 2 higher than other linear or nonlinear functions. The subsequent optimization procedure obtained the variability laws of the parameters a1, a2, and a3 of the seconddegree polynomial functions on the CLWy-GN plane.
Here is the empirical-mechanistic friction model, Equation (8): where, GNnew is the GN value for the new pavement; GNrestart = −0.0121CLWy + 0.905 is the GN value after derubberization; ∆GN' is the loss of non-recoverable friction calculated as the difference between GNnew and GNrestart; ∆GN'' is the loss of recoverable friction calculated as shown in Equation (9): where, a1 = −0.41905GNrestart 2 + 0.6179GNrestart − 0.2299; a2 = −82.445GNrestart 3 + 196.49GNrestart 2 − 155.64GNrestart + 41.03; a3 = −4034GNrestart 4 + 13281GNrestart 3 − 16351GNrestart 2 + 8923.8GNrestart − 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 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 2 higher than other linear or nonlinear functions. The subsequent optimization procedure obtained the variability laws of the parameters a 1 , a 2 , and a 3 of the second-degree polynomial functions on the CLW y -GN plane.
Here is the empirical-mechanistic friction model, Equation (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): ∆G N = a 1 CLW 2 y + a 2 CLW y + a 3 (9) where, a 1 = −0. 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.
the initial post-derubberizing state. Different GN values can be obtained for the same CLWy 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 GNrestart line. In other words, it depends on the history of the accumulated loads and maintenance operations performed.
A graduated abacus on the plane CLWy-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.

Validation Phase
Validation of the empirical-mechanistic model was carried out by comparing the results predicted by Equation (8) (7), see Table 4.

Validation Phase
Validation of the empirical-mechanistic model was carried out by comparing the results predicted by Equation (8) (7), see Table 4. 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.

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. 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 GNpredicted-GNobserved 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.

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. 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 GNrestart value, as shown in Figure 8. In the same way, using the equation shown in Figure 9, it is possible to schedule further maintenance operations.  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 GNrestart 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: In the same way, using the equation shown in Figure 9, it is possible to schedule further maintenance operations.  In the same way, using the equation shown in Figure 9, it is possible to schedule further maintenance operations.

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.