{'fraction_non_alphanumeric': 0.02413047680978576, 'fraction_numerical': 0.011054675829100687, 'mean_word_length': 4.591588050314465, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An efficient physical-based method for predicting the long-term evolution of vertical railway track geometries

The dynamic wheel-rail contact forces resulting from the interaction between vehicle and track are responsible for the local track settlement. If these local settlements vary along the track, geometric irregularities develop further amplifying the dynamic loading of the track caused by the interaction between the vehicle and track. In this work, an efficient vehicle-track interaction (VTI) model is presented for predicting the long-term evolution of vertical track settlement during operation. The VTI model has two interacting components – vehicle and track. The vehicle model describes the vertical dynamics of an 8th of a car. The track model considers an elastic rail on discrete (sleeper) supports. Each sleeper location can have its own stiffness, relative height and settlement characteristics. Dependent on the distribution of stiffness and settlement behaviour along the track together with the initial track geometry, each sleeper settles dependent on the number of load cycles (vehicle passes). The track model is initialized with measured vertical track geometry data and static track deflection data at the beginning (day 0) for two types of track sections in the field, a track section where concrete sleepers with Under Sleeper Pads (USP) are used and a track section where only concrete sleepers are used. Using the same settlement model parameters (constant along the track) for the two tracks, the physical-based VTI model can predict the different track geometry quality evolution for both tracks over 350 days. Finally, the VTI model is used to assess the track geometry deterioration when the track/vehicle properties are changed. The prediction strength of the fast VTI model based on the physical understanding can assist in designing and optimizing tracks and in supporting of maintenance activities.


Introduction
Ballasted track represents the most commonly used type of railway tracks. Rails are discretely supported by sleepers which are bedded in a compacted bed of ballast. 1 The dynamic wheel-rail contact forces are transmitted from the rail to the sleeper, from the sleeper to the ballast and finally from the ballast to the subgrade. These dynamic forces are responsible for the deterioration of the track geometry and related settlement phenomena, 2 e.g. by damaging the ballast, 3 sleeper, and plastic deformations in the subgrade. 4 When these local settlements vary along the track, geometric irregularities develop even further, causing an amplification of the dynamic loading of the track due to interaction between vehicle and track. [4][5][6][7] Therefore, a good physical understanding about the development of such track geometry evolution is of great interest, e.g. to assess the running safety and ride comfort of vehicles and to reduce long-term track geometry deterioration by introducing new track components or improving properties of existing track components, therewith reducing maintenance-related costs. 8 Track stiffness is known to be an important parameter influencing the dynamic interaction between vehicle and track. 9 This physical quantity influences the track quality and its long-term deterioration. 10,11 Track stiffness is known to vary from sleeper to sleeper and over longer lengths of track. 9,[11][12][13] Load distributions on the sleepers are uneven due to variability of the track stiffness. 12 Field measurements have shown an accelerated track geometry deterioration on track sections with high variation in the substructure stiffness created by vehicle-track dynamic interaction. 14,15 Track geometry deterioration can be reduced by using Under Sleeper Pads (USP) at the sleeper-ballast interface. 16 Both, the variations in track stiffness and track geometry deterioration have shown to decrease when USP are used. 13,16 Moreover, the peak sleeper load, responsible for the incremental sleeper settlement normally decreases when USP are used due to a better load distribution. 16,17 Local variations in track stiffness and track geometry contribute significantly to the track geometry deterioration and should be considered in computation models. 8,10,12 Several sophisticated models have been developed in the last decades to describe the short-term vehicle-track coupled dynamics accounting for these local variations, see 2,18,19 and reference therein. Predicting the long-term evolution of the vertical track geometry under operation can be done by simulation models. 2,4,5,7,10,12,20,21 Such a simulation model is used by Grossoni et al. 10 to analyse the effects of parameters such as vehicle unsprung mass, vehicle speed, track stiffness and track irregularity on the evolution of track geometry for the nominal track. Using the simulation models, interesting insights have been obtained on transition zones such as rail joint, 7,22 rail weld irregularity, [23][24][25] culvert transition, 6 bridge abutments, 21 railway turnouts, 22,26,27 among others.
Choice of the settlement model and its parameter in the computation models is crucial for accurately representing the track geometry deterioration process. 6 Dahlberg, 28 Sysyn et al. 29 and Abadi et al. 30 provide a summary of the settlement models from the literature. The loading (force) amplitude is not explicitly considered in most of the empirical settlement models. Laboratory experiments on railway ballast have shown that the ballast settlement is a nonlinear function of the loading amplitude. [30][31][32] Hettler settlement model 33 considers a power of 1.6 of the loading amplitude in the settlement relation. During the service, the loading amplitude of the sleeper is expected to change. Stewart, 31 Ford 34 and Mauer 35 present methodologies to account for the changing loading amplitudes for predicting the track settlement in simulation. In their work, phases of constant amplitudes were assumed in between changing loading amplitudes in their simulation methodology and qualitative insights were obtained into the physical phenomena involved for the track settlement. Varandas et al. 6 proposed a settlement model for an arbitrary loading sequence and accounting for the loading history, and presented a validation using a simulation model for track settlement near the culvert. Nielsen and Li 36 use a linear settlement model with respect to vehicle pass, based on the field measurements for track geometry degradation predictions. Dahlberg 37 points out that linear settlement regime occurs anywhere between 50000 and 1 million vehicle passes.
Though significant work has been done using simulation tools for predicting track settlement, there are still open research questions of interest for the railway community, such as: (i) What is the physical explanation for track settlement -why different tracks subjected to the same loading conditions can lead to different settlement behaviour 37 ? (ii) What is the role of track stiffness on the track irregularity development 10,37 ? (iii) How is the distribution of track support properties linked to the track geometry deterioration rates 14 ? (iv) How can the track irregularity and stiffness variation be included in the track model for making realistic forecast of track performance over a longer track lengths 12 ? Understanding the development of such track irregularities or knowing how the growth of such track irregularities can be reduced is of great interest, especially with respect to reducing maintenance effort and related costs. 16,[38][39][40][41] Here, a simple and fast, physical-based vehicletrack interaction (VTI) modelling approach is presented for predicting the evolution of vertical track geometries. The focus is put into representing the realistic track stiffness and track geometry variation along the track in the simulation model. Computationally efficient models are important to allow for reliable predictions for the long-term settlement behavior of long track sections. This paper presents a balanced approach where different simple (but still physical) sub-models are linked to predict the track settlement behavior in a fast manner. The paper is organized as follows. The next section briefs about the track deflection and track geometry measurements in the field for the two types of track sections with and without USP. The models for the dynamic vehicle and track systems are presented in the Modelling approach section. The settlement model utilized and the simulation procedure for predicting the track geometry deterioration process are also described here. The Results section shows the comparison of track geometry quality parameter evolution observed in the field with the VTI model. Model application to study the effects of track/vehicle parameter on the track geometry evolution is also presented here. Finally, conclusions and outlook are drawn in the last section.

Field measurements
In this work, field dataset provided by Swiss Federal Railways (SBB Infrastruktur) for a tangent track. The track consists of two types of sections: track section where concrete sleepers with Under Sleeper Pads (USP) are used and a track section where only concrete sleepers are used. This allows studying the effects of USP on track geometry evolution, building up on the basis of Discrete Element Method (DEM) study in Kumar et al., 38 where micromechanical effects of USP on ballast behaviour was investigated in a box test. The track is in operation since 2003 and the datasets used here are after the tamping. The measurement performed is between location 43 km and 47 km as shown in Figure 1. UIC60 rail type is used and the sleeper spacing is 0:60 m. Stiffness transition between the two track sections occurs at location 44:73 km. In this work, for simplicity, the stiffness transition is not modelled and both track sections are simulated separately. Track section with USP is taken between 43 km to 44:5 km. Track section without USP is taken between 45:5 km to 47 km. The most representative train speed on the track is V ¼ 200 km=h. The track is loaded with approximately 0:1 mega gross tons ðMGTÞ per day with vehicles of axle loads 10:5 � 21 tons range. The heaviest vehicle running on the track is a locomotive. The measured track properties are discretized and interpolated at sleeper location using FANI toolbox -a Matlab based tool developed internally for railway specific signal processing.

Static track deflection measurements
Figure 1(a) shows the vertical static track deflection w i for each sleeper along the track. w i measurement was done by the track recording car one time in the beginning with M TD ¼ 10 tons 8th vehicle load. The deflection of the track due to load is due to the net rail support modulus consisting of pads, ties, ballast, and subgrade arranged in series. 42 Besides fluctuations, w i is fairly constant throughout each section. Average hw i i for the track sections with and without USP are �1:9 mm and �0:88 mm respectively. Higher track deflection for the track section with USP is expected since the introduction of an elastic layer at sleeper ballast interface reduces the total track stiffness, resulting in a higher hw i i j j. The fluctuations in w i are higher for the track section without USP compared to the track section with USP.

Track geometry measurements
19 track geometry measurements for both left and right rail were carried out over 350 days in the D1 range altogether. The track geometry measurements were done by the track recording car with M TG ¼ 8:25 tons 8th vehicle load using chord measurement technique. 8 For simplicity and due to the design of the quasi-2D VTI model, average track geometry of the left and the right rail is considered as the vertical track geometry in this work. Even though the track geometry that comes from the measurement car is filtered, in this work, this mean-free TG i is considered as the actual mean-free track geometry for the two types of tracks. Figure 1(b) shows the vertical mean-free track geometry TG i at the beginning (day 0) for each sleeper along the track. Compared to the track section with USP, track section without USP has higher fluctuations with higher peaks (maximum/minimum values relative to zero). This first track geometry TG i will be used for track model initialization as discussed in the Model initialization section. Figure 1(c) and (d) show two more mean-free track geometries after 252 and 350 days respectively, for the two track sections. Though the track geometry signals are not synchronized, its evolution can be seen for both track sections, with increasing (maximum/minimum) peaks. This is because the track in operation generally deteriorates with time. Quantification of vertical track geometry quality can be described by standard deviation r and peak-to-peak values. 8,43 Figure 2 shows the evolution of r and peak-to-peak values, measured for the whole track length, for the two track sections. At the beginning (day 0), track section with USP has a better track geometry quality compared to the track section without USP. Besides fluctuations, the track geometry quality falls with the operation (increasing days) for both tracks, as can be seen by the increasing r and peak-to-peak values. Note that the track section without USP shows a faster track geometry deterioration (faster r and peak-to-peak values increase) in 350 days compared to the track section with USP. Figure 3 shows the vehicle-track interaction (VTI) modelling approach followed in this work for predicting the vertical track geometry evolution. The dynamic vehicle and track system is shown in the left part. From the left part, a load history is obtained for each sleeper during a vehicle pass, which is dependent on the distribution of stiffness along the track and track geometry. The load history can be used as an input in the right part of Figure 3 to obtain the incremental settlement behaviour of each sleeper dependent on the number of vehicle passes (load cycles). This is used as an input in the left part of Figure 3 to update the track properties (relative sleeper height) where the next vehicle pass takes place in the next loop. Here, for simplicity, an empirical settlement model is used to represent the right part for calculating incremental settlement. With this modelling approach, the track geometry evolution can be obtained from the VTI model.     Figure 2. Evolution of the track geometry quality parameters: standard deviation r (top row) and peak-to-peak values (bottom row) in 350 days for the two track sections with (left column) and without USP (right column). Fluctuating circles represent field measurements and smooth curves passing through measurement dataset are VTI model prediction results for the full 1:5 km track. The shaded region represents the range of variation when the track geometry quality parameters are obtained from the field measurements in different 200 m subsections and the smooth lines without markers are the corresponding range of variation predicted by the VTI model. . VTI modelling approach used in the present work. Left part represents the dynamic, vertical vehicle and track system. Right part represents sketch for a total sleeper settlement including ballast (taken from Kumar et al., 38 where particle-based DEM model is used to gain deeper physical understanding about the ballast bed behaviour) and subgrade. During a vehicle pass, the force on each sleeper is an output of the left part, which is used as input into the right part to estimate the local track settlement. The local settlement is inserted back into the left part to update the track geometry and the next vehicle pass takes place.

Dynamic vehicle and track model
Both, the primary and secondary suspension systems are considered. k P and c P are stiffness and damping coefficients of the primary suspension system of the vehicle. k S and c S are stiffness and damping coefficients of the secondary suspension system of the vehicle. Gravity g ¼ 9:81 m=s 2 is acting downwards. Vertical motion of the car body is described by: where '_' is the time derivative. Vertical motion of the bogie is described by: and vertical motion of the wheel is described by: where Q dyn is the instantaneous wheel-rail dynamic contact force, determined using the nonlinear Hertzian contact theory for interaction between rail and wheel. A summary of the parameters for the 8th vehicle, representing a locomotive, is provided in Table 1. The considered 8th vehicle has The vehicle moves to the right on the track with representative speed V ¼ 200 km=h.
Track model. The left part of Figure 3 shows also the track model representation used in this work. The rail is simulated as Euler-Bernoulli beam discretely supported at sleeper locations. For the low to medium frequency content considered in this work, it is sufficient to model the rail with the Euler-Bernoulli beam. 5,44 In the beam, each node at the sleeper is connected to the ground using a single springdamper assembly, representing the net stiffness and damping of the track bed below the rail (rail-pads, ballast, sub-ballast and subgrade). Although pads and granular materials present a nonlinear material behaviour, 6,42,45 this effect is not considered in this work for simplicity.
The generalized beam element of the track structure 5 is presented in Figure 4(b). The rail (beam element) is characterized by the cross-sectional area A, second moment of area I, elasticity modulus E and the material density q. EI is the flexural stiffness. The length of the element l is the sleeper spacing. i and j are the sleeper (node) location of the beam element. The momentary downward load Q dyn due to wheelrail contact, is applied at distance x from the left node. k i , k j and c i , c j are the net (linear) spring stiffness and damping respectively at the two nodes. In the presented track model, the height difference (relative height) of the two node locations, i and j, h i and h j can also be accounted. This allows to include the track geometry irregularities in the track model. The simple track model used in this work is inspired by Lei 5 and is modified to account for the track geometry irregularity in the track model.  The nodal displacement and force of the generalized beam element of the track structure are defined as follows: where z i ; z j are vertical node displacements, h i ; h j are rotation angles, V i ; V j are vertical node forces and M i ; M j are bending moments at node i; j respectively. The superscript "T " represents the matrix transpose. Next, different element and force matrices, as well as the link between them are defined.
Element stiffness matrix. The explicit expression for the element stiffness matrix of a generalized Euler-Bernoulli beam element of the track structure is: Assuming that the elastic force is proportional to the node displacement, according to Figure 4(b), the element stiffness matrix induced by the discrete elastic support is: where 1/2 accounts for the sharing of nodes in the current beam element with adjacent beam elements. When assembling the global stiffness matrix, this multiplier vanishes. The net stiffness matrix of the generalized track element is derived by summing up equations (4) and (5): where the subscript "T" represents the track quantity.
Element mass matrix. The enlarged consistent mass matrix of the generalized beam element of the track structure is: Since, the mass of the sleeper and the ballast is not taken into account in this work, the net mass matrix of the generalized track element is: Element damping matrix. The damping matrix of the beam element of the generalized beam element is often expressed as follows: This damping is called proportional damping, with damping coefficients p and q, related to the damping ratio and natural frequency of the system. 25 Just like the element stiffness matrix of track structure, in addition to the proportional damping caused by rail, the damping force caused by discrete supports should also be taken into account, expressed as: The net damping element matrix of the generalized track element is derived by summing up equations (9) and (10): Element force matrix. The load vectors on the track element are: distributed body weight force, vertical concentrated force (wheel-rail contact force), and the stored spring force due to height difference in the nodes of the track element. The element load vector due to distributed body weight force applied on the full length of the beam element is: The element equivalent node load vector induced by vertical concentrated wheel-rail contact force Q dyn acting downwards at a distance x from the left node location, as shown in Figure 4(b), is 46 : Note that f e Q is 0 for all elements where no external vertical loads are acting, i.e. all elements except where Q dyn is acting. The element load vector induced by stored spring force due to height difference h i and h j between the node is: The force element matrix of the generalized beam element of track structure is derived by summing up equations (12) to (14): Dynamic equation for track model. Dynamic equation for track model. The dynamic equation for track model is described as 5 : where are the global mass matrix, damping matrix, stiffness matrix, and global load vector of track structure respectively. a T ¼ z i ; h 1 ; z 2 ; ½ h 2 ; . . . ::; z N�1 ; h N�1 ; z N ; h N � T is the global node displacement vector of the track structure. A summary of the parameters for the track model is provided in Table  2. Note that by modifying and arranging the element mass, damping and stiffness matrix in equations (4) to (15), additional components such as rail-pads, sleeper, ballast can be included easily in the track model.
Interaction model. The vehicle and track systems are coupled through the interaction force Q dyn between the wheel and the rail, determined using the nonlinear Hertzian contact theory: where z W ; z R and g are instantaneous wheel vertical displacement, rail displacement and rail surface irregularity at the wheel-rail contact point respectively. z W � z R þ g ð Þ is basically the indentation of the wheel in the rail surface. The rail surface irregularity g is set to zero in this work. G is the contact deflection coefficient, which for a new wheel with conical thread according to Lei 5 is: where wheel radius R is taken as 0:458 m. The vehicle and track nonlinear system is solved in time domain using the implicit Newmark integration scheme as two separate subsystems. The coupling of the two subsystems through interaction forces is achieved using a cross-iteration algorithm. 5

Track settlement model
A sleeper settlement curve during a tamping cycle is sketched in the right part of Figure 3 consisting of ballast and subgrade settlement. The settlement of ballast during the tamping cycle comprises of three different phases. 38 The settlement is relatively fast in phase I due to the fresh nature of the ballast bed. The settlement of subgrade is generally linear during the tamping cycle. 6 When the service begins normally, the track has already experienced several vehicle passes. Thus, in this work, it is assumed that the available measurements have been done in phase III, where the settlement is slow and linear. 7,21,28 Track settlement highly depends on the size and type of sleepers, quality and the behaviour of the ballast, sub-ballast and subgrade, traffic history, maintenance history, environmental conditions, and the lift given to the track during maintenance. 28,37 Moreover, it is expected that track loading and (peak) force on the sleeper changes continuously during the tamping cycle. Hettler settlement model, 33 with sleeper force as a direct parameter in the formulation, is considered most suitable in this study, described as: where h N i is the sleeper height of the i th sleeper after N vehicle passes, with initial height h 0 i . Thus, h N i � h 0 i is the sleeper settlement after N vehicle passes. F i;max is the maximum force on sleeper in kN. a i ¼ 1:6 describes nonlinear relationship between force and settlement in the Hettler settlement model. 6,30,33 r i ¼ 9:5 � 10 �7 m=kN 1:6 h i and C ¼ 0.43 are empirical constants. 30 In this study, it is assumed that the majority of the track's life-time operation occurs in the linear settlement phase III, the generalized Hettler model in a linearized form (obtained simply as the first derivative of equation (19)) is used, represented as: where Dh i is the incremental settlement (change in h i ) of the i th sleeper during DN number of vehicle passes. In this work, the power exponent a i ¼ 1:6 is used in the settlement model. s i (not known) is a function of ballast and subsoil properties, operational history 28,37 etc. and can be derived from equation (20) from a well-synchronized track geometry evolution measurements from the field. In this work, s i is the only free parameter that is fine-tuned to get a similar track geometry deterioration with days in operation observed in the field, as seen in Figure 2. Both a i and s i values are kept constant throughout the track and in time. The same set of a i and s i values are used for the settlement predictions for the two types of tracks with and without USP. Keeping the same a i and s i for both tracks is a logical assumption, as the track settlement during the (linear) settlement phase III is dominated by wear of ballast stones 38 and linear subgrade settlement, 6 which depends primarily on the subsoil properties, operational history etc. and not much on the elastic layers at the sleeper ballast interface.

Simulation procedure
A flowchart describing the procedure for VTI simulation is presented in Figure 5. The first step is the initialization of the track model. Using the initial condition of the track (vertical static track deflection for a given load and track geometry under loaded conditions), the track model is initialized to obtain the spring stiffness k i and relative height h i for all sleepers along the track (initialization procedure is described in the Model initialization section). The next step in the VTI simulation procedure is running the 8th vehicle with constant speed V on this prepared track. The dynamic vehicle pass is performed using the crossiteration algorithm 5 and the wheel-rail dynamic Q dyn -force is extracted. After the vehicle pass, a time history of the force that each sleeper has experienced during that vehicle pass is extracted as the next step. The peak sleeper force F i;max on the i th sleeper is also extracted during the vehicle pass. The final step of the time simulation loop is calculating the local sleeper settlement due to loading by the vehicle within a given number of DN vehicle passes (MGT). The number of vehicle passes DN in the simulations is set according to the (variable) time duration between track geometry field measurements. F i;max is assumed to be not changing within this time duration for different sleepers to speed up the computation. Using the known track loading MGT per day data, the vehicle load, and time duration between field measurement, the number of vehicle passes DN can be calculated. The incremental sleeper settlement Dh i on the i th sleeper is then calculated using equation (20). A new height h i ! h i � Dh i is given to all sleepers. This is how a new track geometry is obtained after the first vehicle pass in the VTI model on which the next vehicle pass is performed. Repeating this process in the VTI simulation loop gives the evolution of the absolute vertical track geometries. The stiffness change during a settlement cycle is slow for a nominal track 39 and therefore assumed not changing for both tracks during the 350 days. It is important to highlight that the changing stiffness during the tamping cycle can be easily included in the VTI model.

Model initialization
The first step towards the VTI simulation loop is the initialization of the track model to set the local track stiffness and relative height for all sleepers along the track. The model initialization step requires the vertical track stiffness and vertical track geometry field measurements done at the beginning (day 0). During the field measurements, measurement car moves with a certain velocity along the track. However, for simplicity, the dynamic interaction between vehicle and track is assumed negligible in this work, and that the field measurements contain only static track behaviour. The static part of the dynamic track model shown in equation (16) is: where a T and F T are the global displacement vector and global node load vector of the generalized track structure. In other words, the global displacement vector a T can be obtained when the track properties and loading conditions are known by: Local track stiffness: As the first part of the model initialization step, stiffness along the track is extracted. This is done using static vertical deflection measurements w i shown in Figure 1(a). The static vertical deflection w i is independent of the relative height of the sleepers. 47 This is to say that static deflection due to external load at any location on the track is dependent only on the stiffness distribution of the neighbouring sleepers and is independent of the relative sleeper heights. Hence, the global node load vector of the track structure reduces to For the first iteration in stiffness estimation, a track stiffness for each sleeper along the track is assumed using the following relation for Winkler foundation 47 : where w i is static track deflection obtained in the field using measurement wagon of M TD ¼ 10 tons 8th vehicle load. From the 1 k i , global stiffness matrix K T is then evaluated using equation (6). F T is also known at all the sleeper locations when loaded with M TD ¼ 10 tons 8th vehicle load, keeping f e h ¼ 0 for track stiffness measurements. Thus, a vertical node displacement at each sleeper location 1 w i is calculated using equation (23). The obtained 1 w i is compared with the static deflection known from the field w i and a difference norm is calculated. Based on the difference norm, a second stiffness profile is assumed 2 k i along the track using a multivariable optimization toolbox fmincon in Matlab. The process ends when a reasonable tolerance limit in the difference norm k nþ1 w i � w i k < tol is achieved for the static track deflection between estimated value and measurements. Thus, a stiffness profile along the track is obtained from the vertical track deflection measurements.
Relative sleeper height: The next part of model initialization is extracting the relative sleeper height h i along the track using the first vertical track geometry measurements TG i shown in Figure 1(b) (day 0) and the optimized track stiffness k i . TG i is obtained from track recording car with M TG ¼ 8:25 tons 8th vehicle load. In this work, it is assumed that the filtered, mean-free TG i obtained from the field represents the actual mean-free track geometry. TG i is a combination of both: the unloaded track geometry irregularity (mainly due to variations in relative sleeper height h i ) and the deflection of the loaded rail (due to variation in the vertical support stiffness k i ), 47 as also evident from equation (23).
The method of calculating the relative sleeper height hi is straight forward. First an array of N s unknown h i 's with length equal to the number of sleepers N s is created. As evident from equation (23), at each sleeper location, K T (from k i ), f e r , and f e Q (from load M TG ) are known. Still unknowns in equation (23) are a T and f e h (from h i ). Inserting the value of TG i from field measurements in a T for a sleeper location, a single linear equation with still unknown h i 's is obtained. Repeating this process at every sleeper location, N s linear equations are obtained with N s unknown h i 's, which are then solved and to obtain the relative height h i of all the sleepers.
The local track stiffness k i and relative height h i obtained for each sleeper along the track after the initialization step is shown in Figure 6. Figure 6(a) shows the stiffness profile for the two track sections. The average stiffness hk i i in the track section with and without USP are 17:2 kN=mm and 49:8 kN=mm respectively. In other words, hk i i in track section without USP is 2:9 times larger compared to the track section where USP are used. Equation (24) shows a direct relation between average track stiffness and static deflection as k 3 i w 4 i ¼ constant. From this analogy, the average static deflection in the track section with USP should be approximately 2:9 3=4 �2:2 times higher than track section without USP. This is confirmed when looking at the field measurements for w i in Figure 1(a). Higher fluctuations in estimated stiffness k i for the track section without USP is seen in Figure 6(a). This is mainly due to higher fluctuation in w i present in the field for the track section without USP, as seen in Figure 1(a). The relative sleeper height h i obtained after the initialization step shows a variation of �7 mm for both track sections, as shown in Figure 6(b). Track section without USP has higher fluctuations in h i compared to the track section with USP.
To check the accuracy of the obtained track properties k i and h i after the model initialization step, k i and h i , are used in equation (23) to calculate back the static track deflection w i and track geometry TG i . The obtained values are plotted as dots in Figure 1 against the corresponding field measurement dataset. A remarkable agreement with field measurements can be seen for both w i and TG i obtained from the estimated k i and h i . For both tracks, the root mean square error between the field measurements and estimated values for both w i and TG i is less than 10 �4 mm.

Model prediction
The next step in the VTI simulation loop is to make a vehicle pass on the initialized track as shown in Figure 5. Performing VTI simulations for a variety of traffic (axle loads, velocities, etc.) running on the track is computationally expensive. Therefore, a single 8th vehicle, with properties described in Table 1, running at a representative speed V ¼ 200 km=h is used to simplify the computation. All the vehicle traffic variations are somehow fused into the local settlement parameter s i that leads to a similar vertical track geometry quality evolution in simulations as observed in the field.
The 8th vehicle moves from left to the right with speed V ¼ 200 km=h on both track sections with and without USP. Figure 7(a) shows the dynamic Q dyn -force at the wheel-rail interface obtained from the VTI model for both track sections. The wheel-rail contact force has a dynamical behaviour fluctuating around the static vehicle load Q vehicle ¼ M vehicle g ¼ 103 kN. After the vehicle pass, the force-time history and hence the peak force on each sleeper F i;max is also extracted. F i;max for the two track sections is shown in Figure 7(b). The average peak sleeper force hF i;max i in the track sections with and without USP are 33:85 kN and 44:36 kN respectively, also given in Table 3. This is to say that in average, the peak sleeper force in track section without USP is �1:3 times higher than that of the track section with USP. Under a load application, bending of the rail is smaller on a track without USP compared to a track with USP. 17 The force is distributed to fewer sleepers, and thus, higher hF i;max i is expected in a track without USP. Equation (24), though valid for constant track stiffness and static scenarios, can aid to estimate of the ratio of the peak sleeper force between the two types of track sections simulated here: hF i;max i ¼ hk i ihw i i / hk i i 1=4 . The track stiffness for the track section without USP is 2.9 times higher than that of the track section with USP. Thus, the average peak sleeper force hF i;max i in the track section without USP is expected to be 2:9 1=4 �1:3 higher compared to the track section without USP, as confirmed by VTI simulation output in Figure 7(b).
The final stage in the VTI simulation loop is to calculate the incremental sleeper settlement Dh i . Since the track in operation is considered in stage III of the settlement cycle (right part of Figure 3) where settlement occurs linearly slowly, the linear settlement model from equation (20) is used to calculate the incremental settlement. F i;max required in settlement model in equation (20) is taken from the vehicle pass in the VTI model. The number of vehicle passes DN is based on the average MGT track loading and the time duration between the field measurements. Within this time duration F i;max is assumed to be not changing for different sleepers. Thus, the incremental sleeper settlement Dh i along the track is estimated during DN vehicle passes and a new track is generated. On this new track, the vehicle makes the next pass and the same steps as shown in Figure 5 for the simulation procedure of the VTI model are repeated. The repeating steps of vehicle passes, extracting forces and updating track geometry is stopped when a user-defined number of days (or MGT) is reached, which in the present work is 350 days of field measurements.
In the settlement model, a fixed a i ¼ 1:6 for all sleepers is used for both tracks as discussed in the Track settlement model section. The s i parameter is tuned in the VTI model to obtain a similar track geometry deterioration as in field measurements during the 350 days of service. The tuned s i parameter is 1:9 � 10 �12 m=kN 1:6 h i . A relationship between the Hettler to obtain a similar track ers r i and C in equation (19) and the linear settlement model parameter s i in equation (20) can be deduced as: where N � is a characteristic vehicle pass after which the settlement behaviour is linear. Dahlberg 37 points out linear settlement regime after about 50,000 up to 1 million N � characteristic vehicle passes. The value of r i C=N � thus is expected to be 8 � 10 �13 � 4 � 10 �12 ½m=kN 1:6 �, very similar to the s i parameter value tuned in this work. This is an astonishing outcome of the VTI model where the linear settlement model parameter s i is tuned to match the field observation, which turned out to be very similar to the literature data. The peak sleeper force F i;max is higher on the track section without USP, and thus, the absolute track settlement is also higher (data not shown). Figure 8 shows the mean-free, loaded track geometry development TG i for the two track sections obtained from the VTI model for a few operation days as shown in the legend. The same vehicle load M TG of the track geometry measurement car is used to obtain TG i . With increasing time, the VTI model shows lesser TG i development for the track section with USP compared to the track section without USP. Slower TG i development for the track section with USP is mainly due to smaller peak sleeper force F i;max , since all the other parameters in the settlement model are kept the same for the two types of track sections. Effects of initial track geometry and stiffness profile on the evolution of TG i are also important in the development of TG i and will be investigated in future work.
Finally, the vertical track geometry quality parameters standard deviation r and peak-to-peak values are calculated for different TG i obtained from VTI model. Figure 2 shows the evolution of r and peakto-peak values from the field as well as from the VTI model for the two track sections. Since the field measurements at the beginning (day 0) are used for VTI model initialization, r and peak-to-peak values are the same in the VTI model as in the field. r and peak-to-peak values from the VTI model show a nonlinear, plunging track geometry quality behaviour with operation days for both tracks with and without USP. The field measurements show a slower track geometry deterioration (lower slope of r and peakto-peak values curves) for the track section with USP compared to the track section without USP. The VTI model also predicts reasonably well (almost) no growth in the peak-to-peak values for the track section with USP. This is in agreement with the field work from Schneider et al. 16 who also observed a slower track geometry deterioration when using USP in the track. It is important to highlight that the fluctuations present in the field data for r and peak-to-peak values could be because of the measurement uncertainties while the VTI model predicts a smooth track geometry deterioration for both tracks. The fine-tuned s i settlement model parameter could describe well the field observations.
The VTI model was also used to asses the evolution of range of variation for the evolution of the two track quality parameters obtained in different 200 m subsections -a common practice for the infrastructure managers -within the 1.5 km length of the two track sections. Figure 2 shows that the VTI model predicts well the range of variation for r for both tracks, with and without USP. The peak-to-peak value of a small subsection (200 m) is always less than or equal to the peak-to-peak value of a large subsection (1:5 km). The upper range for peak-to-peak value obtained from different subsections is quite close to the peak-to-peak value for the full 1:5 km section as seen in Figure 2. Thus, the VTI model results, based on physical understanding, can describe well the different track geometry deterioration for the two types of track sections which gives further confidence regarding the used modelling approach.

Model application
In this section, the VTI model will be applied to evaluate the track behaviour when the properties of track/ vehicle components are changed. This is particularly important for the vehicle manufacturers and track operators to assess the track geometry deterioration when track/vehicle properties are altered. Two case studies are performed: effects of (i) rail properties and (ii) unsprung vehicle mass on the track geometry deterioration. The idea is to use the VTI model to identify the main mechanisms responsible for different track geometry developments when track/vehicle properties are changed. Using the same initial conditions, static track deflection w i and same initial track geometry TG i as in the field at the beginning (day 0), the track geometry development will be analyzed for the two case studies for both types of track sections with and without USP.
Effects of rail properties. The UIC60 rail is present in the field. Here, the effects of two other types of rail A and B will be investigated with parameters given in Table 3. Only the parameters different to the original UIC60 rail is given in the table and all other parameters are kept the same as of original rail. Rail A has 27% smaller second moment of area (rail bending stiffness) I than the original rail and is similar to UIC54 type rail (slightly smaller crosssectional area A). On the other hand, Rail B has 27% higher I compared to the original rail. The track model initialization step, as done for the original rail in the Model initialization section, is also made for both rails A and B keeping the same initial static track deflection w i and initial track geometry TG i for both tracks. Following the initialization Figure 8. Evolution of mean-free track geometry TG i for a segment in track sections (a) with and (b) without USP after 252 and 350 days, as shown in the legend. Track geometry at the beginning is also plotted for comparison. step, the VTI simulations were performed on these two new tracks made with rails A and B. The same 8th vehicle with properties described in Table 1 runs on these new tracks with speed V ¼ 200 km=h. Figure 9 (top row) shows the wheel-rail dynamic Q dyn -force in the two types of track sections with and without USP for the two rails A and B as well for the original rail. The average Q dyn -force fluctuates around the static vehicle load Q vehicle ¼ 103 kN and no visible difference in the dynamics can be seen for the three types of rails. This has to be expected since the dynamic Q dyn -force is mainly governed by the track geometry TG i variation (in this frequency region) and the vehicle properties, which are kept the same for the three types of rails. A summary of the standard deviation r Q dyn ð Þis given in Table 3 for all cases. Figure 9 (bottom row) shows the peak sleeper force F i;max along the track for the three types of rails and for the two types of track sections. A summary of the average peak sleeper force hF i;max i and its standard deviation r F i;max ð Þ is also given in Table 3 for all cases. Just like in the original case, F i;max is smaller in track section with USP compared to the track section without USP for both rails A and B as well. It can be seen that hF i;max i decreases with increasing I of the rail. The explanation for this behaviour of decreasing hF i;max i with increasing I can be easily deduced from equation (24). For tracks with same initial static track deflection w i ; hF i;max i ¼ k i w i / I �1=3 . Therefore, hF i;max i should decrease with increasing I of the rail, as also observed from the VTI simulations. The standard deviation for peak sleeper force r F i;max ð Þdecreases slightly for increasing I of the rail for both types of track sections. This could possibly be due to the better distribution of dynamic Q dyn -force to more sleepers with increasing I and therefore less fluctuations in peak sleeper force. Figure 10 (top row) shows the vertical track geometry quality parameters r and peak-to-peak values obtained for 350 days service duration for the three types of rails. The VTI model results show a smooth, nonlinear behaviour with the number of days in operation for the different track and rail types. For the same rail type, r and peak-to-peak values increase faster for the track section without USP compared to the track section with USP. This is to say that the track without USP has a faster track geometry deterioration compared to the track with USP. Higher average sleeper forces in the track section without USP compared to track section with USP is the primary reason for higher track geometry deterioration observed for the three types of rail. For both track types, the track geometry deterioration decreases with increasing I value, i.e. rail A with smallest I leads to highest track geometry deterioration while rail B with largest I leads to lowest track geometry deterioration and the original case lies in between these two. This is again linked with the lower average peak sleeper force F i;max observed for higher I. A lower track geometry deterioration with increasing rail bending stiffness I is in agreement with the literature. 28 Effects of unsprung mass. Here, the effects of two other types of vehicles A and B will be investigated with parameters given in Table 4. Only the parameters different to the original vehicle are given in the table and all the other parameters are kept the same as of the original vehicle. Vehicle A has 500 kg less wheel mass M W and 500 kg more bogie mass M B compared to the original vehicle. On the contrary, Vehicle B has 500 kg more wheel mass M W and 500 kg less bogie mass M B compared to the original vehicle. Thus, vehicle A has smaller unsprung mass M W while vehicle B has higher unsprung mass compared to the original vehicle. All the three vehicles have the same vehicle load, which allows investigating solely the effect of unsprung mass on the track geometry evolution using the VTI simulations. Both types of track sections with and without USP from the original  Figure 11 (top row) shows the wheel-rail dynamic Q dyn -force in the two types of track sections with and without USP for the two vehicles A and B as well for the original vehicle. Unlike for the case study for changing rail property I, a clear difference in the dynamic Q dyn -force can be seen here for the three vehicle scenarios considered. For both tracks, the dynamic Q dyn -force fluctuates around the static vehicle load, which is same for the three vehicles. The fluctuations in the dynamic Q dyn -force increases with increasing unsprung mass for both tracks. This is mainly due to the fact that the dynamic Q dyn -force is linked to the track geometry TG i variation (in this frequency range) and the wheel mass M W , which increases with increasing unsprung mass. A summary of the standard deviation r Q dyn ð Þ is given in Table 4 for all cases where it can be seen that the dynamics in Q dyn -force increases with increasing unsprung mass. Figure 11 (bottom row) shows the peak sleeper force F i;max along the track for the three types of vehicles for both tracks. A summary of the average peak sleeper force hF i;max i and its standard deviation r F i;max ð Þis also given in Table 4 for all cases. Just like in the original case, F i;max is smaller in the track section with USP compared to the track section without USP for both vehicles A and B. hF i;max i is quite similar for the three vehicle types for both tracks. On the other hand, r F i;max ð Þ increases with increasing unsprung mass (wheel mass) M W . This is mainly due to the higher dynamics in the Q dyn -force observed for higher M W , as shown in Figure 11 (top row). Figure 10 (bottom row) shows the vertical track geometry quality parameters r and peak-to-peak values obtained for 350 days of service for the three types of vehicles. The VTI simulations show a smooth, nonlinear behaviour with the number of days in operation for the two tracks and vehicle types. For the same vehicle, r and peak-to-peak values increase faster for the track section without USP compared to the track section with USP. Again, this is due to the higher average sleeper force in the track section without USP compared to the track section with USP, observed for the three types of vehicles. For both tracks, the track geometry deterioration is higher when the vehicle with higher unsprung mass M W operates, i.e. vehicle A with smallest M W leads to lowest track geometry deterioration while vehicle B with largest M W leads to highest track geometry deterioration and the original case lies in between these two. This is again linked with the Figure 9. (Top row) Dynamic Q dyn -force at the wheel-rail interface and (bottom row) peak sleeper force F i;max on sleepers along the track during the first vehicle pass for three types of rails, as given in Table 3. Left to right (a thend d-f) is increasing second moment of area I. (Middle column) Original scenario with UIC60 rail, as also shown in Figure 7. Horizontal line is the static vehicle load Q vehicle ¼ 103 kN.  Figure 10. Evolution of the vertical track geometry quality parameters: standard deviation r and peak-to-peak values during 350 days of service for the two types of track sections. (Top row) effect of three types of rails (different second moment of area I) used in the VTI simulations, as presented in Table 3. (Bottom row) effect of three types of vehicles (different wheel mass M W keeping same vehicle load) used in the VTI simulations as presented in Table 4. Fluctuating circles represents the field measurements and smooth curves passing through measurement dataset are VTI model prediction results. Arrows indicate the direction of change in track geometry quality with increasing parameters value. higher fluctuations in the F i;max observed when increasing the unsprung mass, which leads to a higher track geometry development. Higher track geometry deterioration with increasing vehicle unsprung mass is in agreement with the observation from literature. 10

Conclusions and outlook
A new vehicle-track interaction model has been developed for predicting the evolution of the vertical railway track geometry. The model describes the vertical dynamics of an 8th of a car running on discretely supported elastic rails. Each rail support (sleeper) can have its own stiffness, relative height and settlement characteristic. The VTI model initialization is done with the static vertical track deflection and track geometry measurements at the beginning (day 0). After the initialization step, the local stiffness and the relative height for each sleeper along the track is obtained. On this initialized model, the vehicle runs with a constant speed and the wheel-rail dynamic Q dyn -force is extracted. The time-history of the force on the sleeper and the peak sleeper force is also extracted. Each sleeper settles incrementally based on the peak sleeper force, the number of days (vehicle pass/MGT) passed and track characteristics, according to an empirical settlement model in the linear settlement regime. A new (settled) track geometry is obtained and the next vehicle pass is performed. This is how the track geometry evolution is obtained from VTI simulations. The objective of this work is to present a methodology for predicting track geometry evolution. The methodology uses a balanced approach where different physical-based simple sub-models are linked to predict the track settlement behavior in a computationally efficient and reliable manner. The first static vertical track deflection and track geometry measurements obtained in the field were used for the VTI model initialization for two types of track sections: a track section where concrete sleepers with Under Sleeper Pads (USP) are used and a track section where only concrete sleepers are used. Over the 350 days of track geometry measurements in the field, the track section with USP has shown a better global vertical track geometry quality (slower standard deviation and peak-to-peak change) compared to the track section without USP.
On the initialized VTI model, an 8th vehicle with 10:5 tons load runs at a representative speed V ¼ 200 km=h. The incremental track settlement for Figure 11. (Top row) Dynamic Q dyn -force at the wheel-rail interface and (bottom row) peak sleeper force F i;max on sleepers along the track during the first vehicle pass for three types of vehicles, as given in Table 4. Left to right (a the track during the first vehiclunsprung mass (wheel mass M W ). (Middle column) Original scenario with a locomotive vehicle, as also shown in Figure 7. Horizontal line is the static vehicle load Q vehicle ¼ 103 kN. each sleeper is estimated after the vehicle pass using an empirical settlement model in the linear settlement regime. The absolute settlement is found to be higher for the track section without USP due to higher peak forces experienced by the sleepers on this track. Because of higher peak sleeper forces and higher initial geometric irregularities, the overall vertical track geometry quality parameters standard deviation r and peak-to-peak increase faster (higher track geometry deterioration) for the track section without USP. With a fine-tuned settlement model parameter, assuming the same settlement parameters (constant along the track) for both track sections, the VTI model could predict the different track geometry quality behaviours observed in the field for both types of track sections. This gives confidence that the VTI model accounts for the physics in a good way. The VTI model also predicted well the range of variation for the two track quality parameters calculated in different 200 m for both tracks.
Finally, the parameterized and validated VTI model was applied to investigate the track geometry evolution when the properties of track/vehicle components are changed. Two case studies were performed: the effect of changing second moment of area of rail (rail bending stiffness) and the effect of changing vehicle unsprung mass. For both types of track sections, the track geometry deterioration decreases with the increasing second moment of area due to reduction of average peak sleeper forces. On the other hand, the track geometry deterioration increases with increasing unsprung mass (wheel mass), keeping the same static vehicle load, due to increasing dynamics in the wheel-rail contact force.
The first results are very promising and further validation of the vehicle-track interaction model with track measurement data will be done in the next step. The balanced modeling approach presented in this work can answer complex research questions in an efficient way. Thus, the developed methodology can contribute significantly in reducing maintenance activities. Such an approach can be used to optimize the track system, for example in a transition zone between two track sections with different stiffness by defining stiffness properties of elastic layers that reduce the track geometry deterioration. Extending the simple but fast vertical VTI model to lateral direction should be an important addition and will be considered in the future work.