Modelling railway ballasted track settlement in vehicle-track interaction analysis

The geometry of a ballasted railway gradually deteriorates with trafficking, mainly as a result of the plastic settlement of the track-bed (ballast and sub-base). The rate and amount of settlement depend on a number of factors, and for various reasons are difficult to predict or estimate analytically. As a result, various empirical equations for estimating the rate of development of plastic settlement of railway track with train passage have been proposed. A review of these equations shows that they (i) do not reproduce the form of settlement vs number of load cycles relationships usually seen in the field; (ii) do not reflect current knowledge of the behaviour of soil subgrades in cyclic loading; and (iii) are often critically dependent on the curve fitting parameters used, which in turn depend on the circumstances in which the calibration data were obtained. To address these shortcomings, this paper develops a semi-analytical approach, based on the known behaviour of granular materials under cyclic loading, for the calculation of plastic settlements of the trackbed with train passage. The semi-analytical model is then combined with a suitable vehicle-track interaction analysis to calculate rates of development of permanent settlement for different initial trackbed stiffnesses, vehicle types and speeds. The model is shown to be able to reproduce recursive effects, in which a deterioration in track geometry causes an increased variation in dynamic load, which feeds back into a further deterioration in track geometry. The new model represents a significant improvement on current empirical equations, in that it is able to reproduce observed aspects of railway track settlement on the basis of the known behaviour of soils and ballast in cyclic loading.


Introduction
For nearly 200 years, most of the world's railways have run on ballasted track. With trafficking, the geometry of a ballasted track gradually deteriorates, mainly as a result of the plastic settlement of the track-bed (ballast and sub-base). Intuitively, and as demonstrated by experimentation, the rate and amount of settlement would be expected to depend on a number of factors, including [1]:


The loads imposed, including type and amount of traffic (e.g. axle load, speed, vehicle dynamic effects, cumulative tonnage);


The track superstructure characteristics (e.g. rail and sleeper type, sleeper spacing, railpads and any additional resilient layers such as under-sleeper pads), which influence the distribution of loads into the underlying ballast and ground;  The ballast, sub-ballast and sub-grade layer characteristics (e.g. depth, density, stiffness or resilient modulus, ballast specification in terms of particle size distribution and mineralogy, ballast contamination, drainage and pore water pressure conditions), and the ability of these supporting layers to resist cyclic loading.
When geometry defects become too severe, maintenance -usually in the form of automated tamping or manual packing -is carried out to realign the track and enable the continued safe running of trains. Unfortunately, tamping may also disrupt the load-bearing structure [2] and damage individual ballast grains, resulting in a diminishing return period between maintenance interventions until eventually the track-bed requires full renewal [3].
Four major difficulties in predicting the development of settlement are that: 1. On well-performing track, the rate of accumulation of residual (plastic or permanent) settlement with each loading cycle is almost vanishingly small (in the order of a nanometre); classical soil mechanics theories are not well-suited to modelling such small settlements and their gradual accumulation over potentially millions of loading cycles. 2. Settlement may be attributable to either the ballast or the subgrade, and most likely to both [4]. Different types of subsoil and ballast will have different tendencies to settle; even for ballast having the same grading (particle size distribution curve) and mineralogy, the settlement may depend on factors such as the depth of ballast, shoulder slope and the sleeper type and sleeper / ballast interface conditions (e.g. the use or absence of under sleeper pads) [5,6]. 3. Settlement generally arises as a result of both densification (volume change) and lateral spreading (shear deformation) of the ballast and subgrade [5,6]. 4. It is generally differential (rather than uniform) settlements that cause the track geometry to deteriorate to the extent that it needs maintenance [7]. The development of differential settlement can only be replicated in an analysis if there are pre-defined initial differences, for example in trackbed stiffness and / or in loading, that are not usually initially obvious in reality. There is ample evidence from foundation engineering of a correlation between the maximum settlement and the angular distortion -that is, larger settlements generally are likely to be associated with larger differential settlements [8,9].
These difficulties have led to the development of empirical ballast settlement equations as discussed later, although such equations often do not take account explicitly of differences in sub-base, ballast type and geometry, sleeper type or even loading conditions (axle load and speed).
Recent developments in modelling railway track system behaviour have focused on implementing differentially deteriorating track support conditions in vehicle-track dynamic interaction analyses (e.g. [10][11][12][13]), with the short-term dynamics of vehicle-track interaction (e.g. [14,15]) linked to the long-term degradation of the track through an iterative procedure. The approach is usually based on a time domain simulation of vehicle-track interaction in which the force transmitted by the track system (superstructure) to the supporting layers is calculated at each sleeper position and then used as input into a settlement equation for track geometry degradation prediction [1]. However, there is no consensus on which ballast settlement equation to use. This is unsurprising, because the equations are empirical and their applicability depends on prevailing conditions including traffic type, track structure type and ballast condition.
The aims of this paper are to  review current ballast settlement equations, the range of variables and parameters they can take into account, and the conditions for which they were derived  develop an alternative semi-analytical approach to estimating track support system settlement, based on established soil mechanics principles and referenced to field and full scale laboratory test data  demonstrate the implementation of the proposed approach in a vehicle-track dynamic interaction analysis model to calculate rates of differential track settlement and track geometry deterioration, and compare it with previous methods.

Current ballast settlement equations
A prerequisite to an improved ability to predict the development of differential settlement along a section of track is a better understanding of the relationship between plastic settlement and loading, based on the relevant properties of the ballast and the subgrade at a local scale. Ideally, a ballast settlement equation would be able to account for the effects of the 1. number and magnitude of load cycles 2. train speed (allowing for dynamic load) 3. subgrade separately from the ballast 4. condition of the ballast, subgrade and the interface with the track.
Many, if not all, of these factors are acknowledged within track geometry prediction tools used by industry to plan maintenance over route-scale lengths of track. For example, the T-SPA module within VTISM [16] modifies the basic ballast settlement equation through two main factors. These are the "Local Track Section Factor" (LTSF), which scales the empirical track geometry deterioration equation to the locally measured rate; and the "Ballast Condition Factor" (BCF), which attempts to replicate the observed reduction in the time interval between successive maintenance tamps, generally held to be the result of fines from tamping-damaged ballast grains filling the voids. There are also equations in the literature that attempt to link the development of track geometry deterioration (standard deviation from the desired level) to these factors directly (e.g. [17]). However, these are few and are outside the scope of this paper, which focuses on estimating the rate of track settlement at the level of individual sleepers.
Most authors, including Sato [17] and Dahlberg [18], speculate that there are two major stages of ballasted track settlement, occurring as a result of cumulative loading expressed in million gross tonnes (MGT) or, more usually, the number of cycles of an (often assumed constant) load: 1. Stage 1: after tamping, settlement occurs initially relatively rapidly with MGT or number of cycles as the ballast grains rearrange to establish a structure capable of carrying the applied external loads. This is characterised by a reduction in void ratio and a densification of the ballast [17], i.e. volumetric effects dominate 2. Stage 2: settlement occurs at a slower rate, increasing approximately linearly with MGT or number of cycles. This is attributed to a variety of causes, including the lateral movement of the ballast, penetration of the ballast into the subgrade, and ballast grain breakage and abrasion [18]. Essentially, non-volumetric effects (including shear / lateral spreading) dominate.

Basic equations
Most of the equations that have been proposed to characterise track settlement are empirical. Many take one of two forms: 1.
The more commonly used or cited track settlement equations are summarised in Appendix 1; some of these were reviewed in [18]. A review of empirical permanent deformation models for soils in the context of pavement and railway design was recently presented in [27]. Many equations relate the settlement after N load cycles to the settlement after the first cycle, either explicitly or by inference. Arguably, this is a way of taking into account at least implicitly a range of factors including the load per cycle, assuming that this remains constant. Other equations take into account of a number of factors explicitly, such as the stiffness, condition or nature of the ballast and in some cases the subgrade.
Three problems with these types of equations are that 1. they do not reproduce the form of settlement vs number of load cycles relationships usually seen in large scale laboratory tests or in the field 2. the outcomes are critically dependent on the curve fitting parameters used, which in turn depend on the circumstances in which the calibration data were obtained 3. those expressed in terms of the number of loading cycles do not generally account for the effect of differing axle loads. (Equations expressed in terms of MGT may do, but assumeprobably unrealistically -that the effect of an increase in load is linear).
The second of these points is illustrated in Figure 1, which compares the settlement calculated after 100,000 cycles using the example ballast settlement equations indicated in Table 1 with the curve fitting parameters proposed by the original authors, categorised according to the type of test on which they are based.   Figure 1 shows large discrepancies between the calculated settlements, both within and between each category of experimental basis. In most cases, the discrepancies between categories are intuitively unsurprising. The in situ settlements are largest, but potentially include a contribution from the subgrade, which none of the laboratory measurements do. The equations based on triaxial tests simulating a 20 tonne axle load and the full-scale ballast tests are reasonably consistent; both include only the ballast settlement and the equivalent axle loads are the same. The equations based on parameters from triaxial tests simulating a 32 tonne axle load give greater settlements than those based on triaxial tests simulating a 20 tonne axle load, which again seems intuitively reasonable. The only counter-intuitive difference between categories of equation is that the calculations based on reduced scale ballast box tests seem rather high.
Within an individual category, the variation is greatest for the equations based on triaxial tests simulating a 20 tonne axle load and full scale ballast box tests. It is not clear why this should be, but the potential variability of test specimens, the number of datasets involved and factors such as the frequency of loading could all potentially have an influence. Figure 1 highlights the important influence of the conditions under which a settlement equation and its associated parameters have been derived. Harmonisation would require the development of either a common test procedure that takes into account the effects of vehicle-track interaction, or a modelling approach built up by considering the fundamental behaviour of each of the system components.

More complex equations
In tests carried out to millions of cycles, data of track settlement vs number of load cycles generally show two inflexion points [6,35,36] as indicated in Figure 2, which also shows approximately the phases identified by Sato and Dahlberg. A simple logarithmic or exponential relationship is unable to reproduce this form of curve. (Interestingly, data from a large scale test rig simulating 4 million cycles of high speed train loading presented by Zhang et al, 2019 [37], do not show the second inflexion point and do conform reasonably well to a simple Shenton-type logarithmic equation).
To improve the fit with experimental data exhibiting both phases of behaviour characterised by two inflexion points on the logarithmic graph, more complex equations are needed. Figure 2 shows curve fits from two empirically determined equations able to match a second downturn in the settlement rate on a log scale. The first of these (referred to as the Okabe equation or fit) was originally proposed in [38] (in the same form as that more recently presented in [17]); it combines logarithmic and linear terms (Equation 1): where SN is the settlement at cycle N, N is the number of loading cycles (or cumulative MGT), and C1, C2, α and β are empirically-determined coefficients.
A more complex equation was proposed in [39], based on vibration experiments on columns of confined glass particles, and was fitted by [40] to ballast settlements measured in ballast / sleeper tests carried out to a relatively low number of cycles (10,000), in which the ballast was confined horizontally (Equation 2).
where S1 is the settlement measured after the first cycle and S2, , N0 are curve-fitting parameters with arguable physical meanings.
Equation 2 was developed for settlement occurring solely as a result of densification and is not able to reproduce the later downward phase shown in tests to millions of cycles, which is probably attributable mainly to lateral ballast spreading at the sleeper ends / ballast shoulder. To allow for lateral spreading at higher cycles, Equation 2 can be modified by adding the βN term from Equation 1 to create Equation 3: this is the second equation plotted in Figure 2, referred to as the modified fit.
Using algorithms to determine the constants, Equation 3 can be fitted to data from laboratory sleeper settlement tests to 3 million cycles ( Figure 2). However, the large number of constants of Equation 3 means that there is no unique best fit solution, so it is necessary to constrain the permissible ranges of some of the constants. Different stages of settlement are apparent in both the data and the fitted curves shown in Figure 2, with Equation 3 achieving very good fits. However, whether the components / constants of Equations 1 and 3 individually represent different stages of ballast settlement, and indeed the underlying mechanisms responsible, remains a matter of conjecture. A further feature of these equations (and most if not all of the others reported in Appendix 1) is that they are not inherently dimensionally consistent; hence the curve fitting constants have to have units (dimensions) that make them so.

Equations suitable for use in VTI modelling
Vehicle-track interaction (VTI) modelling offers an efficient way of calculating dynamic wheel-rail contact forces and the displacement, velocity and acceleration of each part of the vehicle and the track system as trains travel along the track. Usually, the vehicle is described as a multi-body system with lumped masses, linked through linear or non-linear force connectors (usually, springs and dashpot dampers). The track is modelled using finite element (FE) theory for the rail as a beam and concentrated masses for the sleepers, with the subgrade support represented by springs and dashpots. Different support system properties may be specified for each sleeper to represent longitudinal variation of support stiffness [41]. Short to medium (1 m to 70 m) wavelength rail roughness, relevant for VTI analysis at train speeds up to about 200 km/hr, is modelled on the basis of actual track recording car measurements along a particular stretch of railway. Different initial horizontal levels may therefore be specified at each sleeper position. In the calculation of long term settlements using VTI analysis, each sleeper position can be adjusted as the track irregularity grows, to capture the effects on vehicle dynamics and wheel-rail interaction forces. It is usually computationally too expensive to change the levels of the sleepers following every train passage; developing differential settlement may be replicated adequately by changing the sleeper levels, for example every 1000 train passes.
Settlement equations suitable for use in VTI modelling should be able to reproduce the main observed features of behaviour, without having so many curve-fitting parameters as to result in a non-unique fit. A further requirement is that the input parameter(s) for the settlement equation, such as load applied and track system response (displacements, accelerations, forces, etc.), can be calculated by the VTI model at each loading cycle, capturing both the traffic characteristics and the longitudinal variability of the track. Not all of the equations in the literature fit these criteria; many incorporate these factors within fixed parameter that have been calibrated for a particular case Phase 1 Phase 2 study. The equations proposed by Guérin [28], Sato [42] and Fröhling [34] have been selected for further study as they do meet the criteria. They are described briefly below.

Guérin's equation
Guérin [28] carried out an extensive series of tests at a loading frequency of 30 Hz, representing approximately the tenth harmonic of the car passing frequency at a train speed of 215 km/h. (For a discussion of the relationship between frequencies of loading, train speed and train geometry, see [43] or [44]). The results suggested that the rate of accumulation of permanent settlement of the ballast may be expressed as a function of the current maximum sleeper deflection: where N is the number of cycles and db,max the maximum elastic sleeper deflection. db,max is a function of N, and expression of Equation 3 in non-differential form would require the ability to specify and integrate this function.
The numerical values of the parameters in Equation 3 proposed by Guérin were for the particular circumstances of the tests used to obtain the baseline data. More recently an extension in the applicability of the formula to train speeds in excess of 350 km/h and "normal" or "soft" soil has been proposed [45]. A major drawback of the Guérin formulation is that it is ill-adapted to a change in track configuration, because the ratio of plastic deformation to maximum displacement will likely change with varying ballast type, depth, etc. There is also a dimensional inconsistency owing to the exponent applied to db,max. .

Sato's equation
Sato [42] presented a two part equation (Equation 5 & 6) for settlement as a function of the number of load cycles N and the ballast pressure P or ballast force F. The settlement SN depends on whether a threshold stress Pth has been exceeded: where the coefficients a and c depend on the power b. The coefficients d and e and the pressure threshold Pth depend on the ballast layer thickness (Table 2 and Figure 3).  Sato's use of different equations above and below the threshold strain reflect the importance of stress dependent non-linear soil behaviour. However, the equations contain dimensional inconsistences and the settlements calculated are sensitive to the stress threshold used. The threshold stress must be assumed to increase with stiffness, otherwise a stiffer trackbed support leads to higher stress transfer from the track superstructure onto the ballast and consequently higher calculated settlements, which is contrary to general experience [46]. The main input parameter for Equation 4, the stress on the ballast, may be computed at each step in a vehicle-track interaction analysis.

Fröhling's equation
By analysing measured data of average track settlement with accumulating traffic, Fröhling [34] developed an expression that, while logarithmic in form [10], also takes account of the deviatoric stress at the sleeper-ballast interface by adjusting the local stiffness, and dynamic load amplification at the wheel-rail interface (Equation 7): where SN is the total settlement after N loading cycles;  The input variables Pdyn and k2mi can be determined as outputs from each step of a VTI analysis, but Equation 7 does have dimensional inconsistencies.

Proposed semi-analytical approach
The review of equations currently used to estimate railway track settlement are mainly empirical; do not incorporate a subgrade behaviour consistent with modern soil mechanics understanding and principles (a point also evident in [27]); do not generally account for differences in axle load; and do not reproduce patterns of settlement with trafficking seen in the field. In this section, a new model for calculating ballast settlements is proposed, which is based on the observed behaviour of soils in cyclic loading, is sensitive to axle load as well as cumulative loading, and is able to reproduce the observed behaviour of ballast under cyclic loading.
In the following, is the vertical stress and the total strain = + comprises an elastic (reversible) component and a plastic (irreversible) one . Using a superposed dot to denote increments of a quantity and assuming linear elasticity, elastic strain increments are calculated as ̇=̇, where is the elastic modulus. We define the plastic strain increment ̇ directly, using a rate equation that fulfils the following requirements:  ̇ is an increasing function of , so that larger stress increments cause larger plastic strain increments  Failure occurs at a limiting or ultimate value of stress ; we assume ̇→ ∞ as → .  Plastic strain develops only during loading, and only if the stress exceeds a threshold ≤  ̇ is an increasing function of the difference ( − ), so that the greater the margin by which the threshold stress is a exceeded, the larger the corresponding plastic strain increment.
An expression satisfying these requirements is: where is a material parameter with units of stress, which is interpreted as a plastic modulus. This expression was arrived at heuristically, and is potentially the simplest one satisfying the above requirements while introducing only one additional parameter.
Consecutive load cycles of equal magnitude are known to result in progressively smaller increments of plastic strain, hence must also increase with loading. It is therefore assumed that is a function of an internal parameter , loosely quantifying how the properties of the material change during the course of deformation; it is assumed that ≡ as a first approximation. However, ≤ always; hence must increase at decreasing rate with . In addition, other things being equal, a higher elastic modulus would be expected to be associated with a higher and a higher .
A general expression that satisfies the above requirements is: This is reasonable, in that a higher peak strength and stiffness are both characteristics of dense or overconsolidated soils. In other words, the stiffness and the strength of a soil would be expected to be correlated, with a stiffer soil also being stronger. For any given soil, the value of the C parameter could be determined on the basis of triaxial tests at different void ratios or preconsolidation stresses, as appropriate for the type of soil and the relevant field conditions. Informed by Equation 9 for , the two equations for ̇ and ̇ can be integrated numerically using the trapezium rule to determine the elastic, plastic and total strain response corresponding to any given stress history. Integrating over a typical load cycle 1 → 2 → 1 using 20 increments was found to provide good resolution and accuracy for the stress-strain response.
Although integrating each load cycle is not onerous, in the sense that results for tens of thousands of cycles on a single sleeper can be produced within a few minutes on a desktop computer, implementing this as part of a vehicle-track interaction analysis over a length of track with hundreds of sleepers leads to a significant computational burden. It is desirable to carry out the integration more efficiently, if possible within a single computational step for each load cycle, even at the expense of a loss of accuracy. Assuming that a load cycle 1 → 2 → 1 increases the plastic strain from 1 to 2 , for known 1 , 2 can be calculated as:

(Equation 13)
Using this approximation results in a loss of accuracy in the order of 5%, which is more than compensated for by the observed two orders of magnitude increase in the speed of calculation. Figure 5 shows, as an example, the form of the cyclic stress-strain behaviour calculated using provides an approximate solution, capture the typical behaviour of a soil-type material in cyclic uniaxial loading. The detail of the behaviour might vary at high frequencies and displacements if inertial effects (accelerations) become significant, or evolve with load cycling as a result of material degradation e.g. grain breakage. These issues have not been explored in this paper (and indeed are generally neglected in the literature), but could in principle be taken into account to some extent through suitable testing and parameter value selection.

Implementation of settlement equations in vehicle-track interaction analyses Overview
The methodology adopted to link the short-term (during train passage) and long-term (permanent settlement) ballast behaviour through vehicle-track interaction analysis was described, using the Sato and Guérin equations, in [12], and is summarised in Figure 6. Here, it is expanded to encompass the proposed semi-analytical approach and the Fröhling equation. Results using all four approaches are then compared using data from a case study with varying support stiffness [12]. The vehicle-track interaction model is based on a finite element (FE) description of a rail as a beam on discrete supports at each sleeper location using support characteristics (linear stiffness and damping terms) that can vary along the length of the track. A quarter vehicle, with half a bogie and two wheels, is considered. Wheel-rail contact is described using non-linear Hertzian theory [41]. The numerical model is integrated in the time domain to calculate train movement and track reaction force.
The dynamic contact loads (hence pressures for the semi-analytical equation) between each sleeper and the ballast (for Sato's equation), the wheel-rail contact force (for Fröhling's equation) and the maximum sleeper displacements (for Guérin's equation) were determined in the time domain for the entire duration of the vehicle passage, and then used as input to the relevant settlement calculation. The track level is quantified at each sleeper position from the incremental settlement ΔS for the traffic loading step ΔN. The vertical rail geometry is updated to include the incremental settlement at each sleeper location and the next dynamic simulation is executed, until the maximum number of cycles is reached. The wheel-rail force and sleeper-ballast interface forces then follow from the multi-body dynamics VTI calculation for the given track geometry configuration (rail irregularity plus sleeper settlement). The parameter ΔN, which defines how many simulations are carried out for a given total number of cycles N, is chosen on a case-by-case basis as a trade-off between accuracy and computational costs. The total settlement is evaluated as the sum of the settlement reached at the previous iteration and the incremental settlement in the current step of N cycles.
Parameter values for the semi-analytical model were determined with reference to single sleeper settlement tests in the Southampton Railway Testing Facility (SRTF) [6]. For the other empirical models, the parameter values specified in the original references were adopted. For the avoidance of doubt, the equations or methods investigated and the parameter values applied in this comparison work are summarised in Table 3.

Simplifications
Distribution of the load by the sleepers into the trackbed depends on the applied load, the bending stiffness of the sleeper and the rails and the effective support stiffness, which may vary along the sleeper length. Stresses tend to concentrate directly under the rail seats [20], but the effective sleeper-ballast contact area depends on a number of factors including the ballast specification, sleeper material and the presence of additional layers such as under sleeper pads [47]. To avoid the difficulties associated with defining the exact nature of the sleeper-ballast interface, the stress on the ballast, the deflection and the load are considered herein as either averages over the sleeper footprint or as the gross value applied to the sleeper. Using the VTI model, it is possible to calculate the force at the rail-pad and the force at the ballast surface at each support position considered. The stress averaged over the total sleeper / ballast contact area (per half sleeper or sleeper end) is then calculated as: where Fb is the force (from a single rail) at the ballast level and As is the area of the sleeper soffit.

Results from two case study sites
Case study 1 "Site B" in [12] has been used in this case study to show the influence of the trackbed stiffness on the long-term ballast behaviour. The original dataset was modified slightly to give stiffness values within the range 60 to 132 MN/m, consistent with Fröhling's assumptions (Figure 7). Two half vehicles typical of the site, i.e. a Class 91 electric locomotive and a laden freight vehicle, were modelled using the parameters given in Table 4.  In total, eight simulations (2 vehicle types  4 settlement models) were carried out. 2 For the original dataset, the mean trackbed stiffness was 110.4 kN/mm (0.6% difference) Figure 8 shows, for the Class 91 vehicle, an example of the evolution with increasing number of load cycles of the calculated track settlement below each sleeper (Figure 8 a & c) and the resulting distribution of ballast forces (Figure 8 b & d), according to the Fröhling (Figure 8 a & b) and the semianalytical (Figure 8 c & d) settlement models. In this case the semi-analytical approach gives much faster settlement rates, with peak values after 320,000 cycles being three (mean) to six (maximum) times those of Fröhling. In the modelling reported in this paper, the increased settlement is fed back directly as an increased track roughness at rail level; this is the approach most commonly adopted in the literature -see for example [34,48,49]. However, it leads to much greater dynamic load amplification and variation in load from sleeper to sleeper as indicated by the increased noisiness in forces seen in Figure 8 d for the semi-analytical model. In contrast, Figure 8 b (using the Fröhling equation) remains practically unchanged through increasing settlement cycles. The increased variation in dynamic forces feeds back into further increased differential settlement and localised defects, indicated by the series of deep troughs that initiate and growth with trafficking. Although the calculated track irregularity remains well within current maintenance intervention limits, even for Category 1 track in the UK, it probably highlights a shortcoming of the approach currently adopted in most coupled VTI settlement models of translating the settlement values as rail level irregularities, rather than introducing them as non-linearities between the sleeper and the ballast surface as in [7].
(e) (f) Figure 8 (e & f) shows that both approaches calculate lower settlement rates at locations with higher trackbed stiffness and vice versa. This is in line with expectations, as the rate of settlement calculated using Fröhling's equation is inversely proportional to the trackbed stiffness, while in the semi-analytical model the threshold stress is proportional to the trackbed stiffness. Settlements calculated using the semi-analytical approach show increasing scatter between individual sleepers with increasing number of cycles, owing to the amplification effect of such track irregularities on the dynamic loads as discussed above.
(a) (b) Figure 9 compares the development of the mean settlement with number of load cycles calculated by all four of the models, for both vehicle types considered. There is no significant difference between the settlements calculated using the Guérin and the Fröhling equations for the two vehicle types, i.e. these methods do not seem to be sensitive to the type and speed of vehicle and its dynamic interaction with the track. Furthermore, the settlements calculated using both of these approaches are very small -less than about 0.2 mm using Fröhling and less than 10% of this using Guérin, after 500 000 cycles. Sato's equation and the semi-analytical approach do show a difference between the responses to the two vehicle types. The higher speed and slightly lower axle load associated with the Class 91 locomotive appear to induce greater settlements than the heavier but slower freight wagon, probably owing to the load amplification discussed previously being greater at higher speed. In terms of magnitude, the calculated settlements of up to about 2 mm after 500 000 cycles are in the lowest quintile of the calculation methods based on full scale or field data shown in Figure 1, and about half those in the rig test data in Figure 2 for a 20 tonne axle load. These are at the low end of the generally reported range, but not unreaslistically low as is the case for the Guérin and the Fröhling equations.
In the case of the semi-analytical model, the load amplification calculated for the Class 91 locomotive is sufficient to cause the specified ultimate stress to be reached, and the relatively large increase in settlement with number of loading cycles is associated with failure of the trackbed. This will need to be refined in future studies to ensure that failure point is not reached unrealistically under expected track and vehicle running conditions.

Case study 2
Case Study 2 is based on the site presented in [50], with the vehicles indicated in Table 4 and the trackbed stiffness characteristics summarised in Table 5 and Figure 10. The trackbed stiffness distributions at the two case study sites are compared in Figure 11. The track stiffnesses at the second case study site are generally less than 60 MN/m/se, outside the range of applicability of Fröhling's equation, which is therefore not considered in this simulation.  The trackbed at Case study site 1 is on average about three times stiffer than at Case study site 2, with less variability. These differences are likely to be at least in part a consequence of the different trackbed stiffness measurement and calculation techniques used, as well as differing site conditions. At site 1, falling weight deflectomer measurements were made on unclipped sleepers, whereas at case study site 2 the effective trackbed stiffness was evaluated from analysis of sleeper movements during train passage and could have included the influence of how well individual sleepers were supported while connected to the rails. The FWD method is arguably an evaluation of the best case trackbed support whereas the latter method with trains present is a more accurate representation of how the track support system responds during train passage. While important, evaluating the most suitable measurement technique is outside the scope of this paper. Case study site 2, with a lower modelled trackbed stiffness, yields a higher mean settlement with greater variation (higher standard deviation) than Case study site 1.
(a) (b) Figure 12 shows the curves of average settlement vs number of load cycles for both vehicle types, calculated for the stiffness distribution of Case Study 2 using the Guérin, Sato and semi-analytical approaches. In this case, all three of the approaches considered show some differences between the calculated settlement curves for the two vehicle types, with the slightly heavier (but slower) freight vehicle giving rise to greater settlements than the Class 91 locomotive. The difference is particularly pronounced for the semi-analytical approach, which gives the smallest settlements for the Class 91 locomotive and the largest (up to the point of effective failure) for the freight wagon. It may be that the contrast between the modelled behaviour of the two sites, with speed being apparently more damaging than load at the first and vice versa at the second, is a result of the different trackbed stiffnesses (the first site being on average three times stiffer than the second). As in the first case study, the displacements are generally at or below the low end of the expected range. This illustrates that further work is needed to understand the physical significance, and obtain representative values of, the parameters used in the equations and models.
The semi-analytical model would benefit from being refined through the selection of parameters more suitable to each case study site. Nonetheless, it does have the potential to capture a range of different track behaviours that the more empirically-based equations do not. Realistic representation of the track in vehicle-track interaction analysis, and the calculation of track cumulative settlements as a result of train passage is in its infancy. The purpose of this paper has been to demonstrate the feasibility and potential of a soil mechanics based approach. Further work is needed, both to understand the physical meaning and quantify the parameters underlying the new model; and to validate the approach with reference to high quality, long-term datasets of track settlement, which at present are rare to non-existent.

Summary and conclusions
Equations that have been proposed to represent the gradual accumulation of plastic settlement of railway track with train passage have been reviewed. Three were selected for further study, on the basis that they are able to reproduce the main observed features of track settlement behaviour, and that they can take as inputs the outputs from an associated vehicle-track interaction model.
A semi-analytical expression, based on the known behaviour of granular materials under cyclic loading, has been developed. This expression is able to reproduce the accumulation of plastic settlement with each load cycle, with the amount of plastic settlement per cycle related to the stress in excess of a threshold stress. The threshold stress increases with the number of load cycles (work hardening), and with the initial stiffness of the trackbed. It also features an ultimate stress, at which plastic deformation continues unchecked. When combined with a suitable vehicle-track interaction analysis, the semi-analytical model has been shown to be able to capture differences in the rate of development of permanent settlement as a result of differences in the initial trackbed stiffness, vehicle type and speed; and the development of rail roughness through differential settlement. It is also able to reproduce recursive effects, in which a deterioration in track geometry causes an increased variation in dynamic load, which feeds back into a further deterioration in track geometry.
While there remains considerable scope for refinement of the semi-analytical approach, particularly through the selection of parameter values appropriate to different site conditions, it shows considerable promise in being able to reproduce observed aspects of railway track settlement on the basis of the known behaviour of geomaterials (soils and ballast) in cyclic loading.
The work has also highlighted areas in which the vehicle-track interaction modelling approach needs improvement, in particular by applying the settlement growth at the interface between sleepers and the ballast, rather than as an equivalent rail irregularity as adopted in much of the current literature. Evaluation of the trackbed stiffness for input into such models also needs further research.