Geomechanical Modelling of Railroad Ballast: A Review

Traditional ballasted tracks have been used intensively around the world with ballast as the main material for tracks. Ballast has a significant contribution to the track alignment, stability and sustainability. After service, ballast deforms and degrades. Periodic ballast maintenance is needed which is a time and cost expensive activity. Understanding the mechanical behaviour of railroad ballast leads to better design and efficient maintenance. From the literature, there are two main approaches used to understand the mechanical behaviour of railroad ballast; large scale experimental and modelling. This paper aims to review the state of the art of literature on the modelling approaches used to understand ballast mechanical behaviour. It discusses the key findings from each modelling approach in understanding ballast mechanical behavior. It presents the main concerns and limitations of each modelling approach from different perspectives related to ballast modelling. It summarizes the limitations, gaps and gaps’ developments of the researches used to understand ballast behaviour via modelling approach.


Introduction
The main role of a railway is to provide a mean to transfer people and/or goods between stations in a safe and sustainable way. A railway track is considered as the basic element of railway infrastructure. It provides a supportive platform to transfer train loads to the ground. There are two main systems of railway tracks; ballasted and ballastless systems [1]. The low cost and higher experience of ballasted tracks are the main reasons of the large number of ballasted tracks worldwide with relative to ballastless tracks [2]. Ballasted tracks have been introduced first to railways while ballastless tracks have been introduced in the 1960s [3].
The two main structures of ballasted track are superstructure and substructure [4] as shown in Fig. 1. The basic element of ballasted track is the ballast layer. Many researchers [4][5][6][7] emphasized on the significance of ballast layer and summarized its main functions as follows: transmits the loads from superstructure to the ground and acts as a bearing platform; enhances the track stability; reduces noise and vibration by providing resilience and damping; provides the track a good level of voids to lose dirt; resists the frost action; reduces vegetation.
After service, ballast material loses its functionality. Ballast deforms and degrades. Ballast maintenance is required which is an expensive activity [8]. The importance of the ballast layer and costly maintenance raise the interest of researches about understanding the mechanical behaviour of railroad ballast. Understanding the mechanical behaviour of railroad ballast results in a better ballast design and efficient maintenance. From the literature, there are different approaches used to understand railroad ballast behaviour; field tests, large scale experimental tests and modelling approach.
Field tests are used rarely to understand the mechanical behaviour of railroad ballast [9][10][11][12][13][14]. Field tests are risky, difficult to control, difficult to collect data during operational time and may require traffic control or stop [15].
There are number of traditional lab tests used to identify the mechanical and physical characteristics of granular material like conventional triaxial test, conventional direct shear test, petrographic analysis, crushing test and Los Angles abrasion. However, Indraratna et al. [8] recommend avoiding the use of conventional tests for understanding the mechanical behaviour of granular material as they rottenly produce confusing results due to the large granular particles size with respect to test sample size. McDowell et al. [16] showed the need of large scale equipment to experiment ballast strength. That raise the requirement of the lab 1 3 experimental large-scale tests to be used in order to understand ballast mechanical behaviour.
Alternatively, modelling the mechanical behaviour based on theoretical models that reflect the real mechanical behaviour of railroad ballast is introduced. The use of modelling approach in understanding the mechanical behaviour has been increasing over the years, especially with the enhancements of computational resources. Modelling approaches can understand the behaviour of railroad ballast under different parameters related to loading conditions and material properties. Modelling approaches develop the knowledge and tools needed for predictions. Such tools have the potential to make huge financial savings within the design constrains, besides, the good understanding and knowledge of the material beahviour under different conditions.
There are two different types of modelling approaches used in the literature to model the mechanical behaviour of railroad ballast based on the problem-solving type; analytical and numerical.
This review paper aims to review the literature on understanding the mechanical behaviour of railroad ballast using modelling approach. The paper discusses the advantages and disadvantages of different modelling approaches, besides, the gained knowledge of ballast behaviour understanding (macroscopic and/or microscopic) from each modelling approach. It highlights the different standpoints of ballast modelling like particle shape, particle size, contact detection models, material models and software packages. Moreover, it summarizes the limitations, gaps and gaps' developments of the researches related to ballast material modeling.

Analytical Modelling
Analytical models of ballasted tracks are founded on mathematical models. Those models describe the mechanical behaviour of each track element based on set of equations related to equation of motion.

Origin and History
In 1867, modelling of ballasted track was addressed by Winkler [56]. Which is known as Beam on Elastic Foundation (BOEF). In 1946, Hetényi [57] illustrated Winkler's model that Winkler modeled the ballasted track as a single beam on elastic homogeneous foundation. Winkler beam has been modelled for different beam theories. Mathews [58,59] modeled Euler-Bernoulli beam on infinite elastic foundation with and without viscous damping. Chonan [60] used Timoshenko beam theory to model the beam on elastic undamped foundation. Kim and Roesset [61] used structural damping in his model to take into consideration the energy losses. Bogacz et al. [62] compared between the responses of Euler-Bernoulli beam and Timoshenko beam models.

Ballast Element Model
From the literature, the ballasted track elements have been analytically modelled using different structural elements based on their characteristics [63]. In general, ballast layer is modelled as mass-spring system with viscous damper to count for the damping effect or with structural damper to count for the energy losses. Sleepers are modelled as rigid  [4] masses. Rails are modelled as beams based on two beam theories; Euler-Bernoulli and Timoshenko beam theories [64,65]. Track subgrade is modelled either rigid foundation or half space.

Analytical Modelling Methods
There are two methods to analytically model ballasted track based on support nature; continuous or continuous discretely supported models as shown in Fig. 2.
In the first method, rails are modelled as beams on a continuously supported foundation. Although, the first approach does not show all the features of ballasted track such as sleeper spacing, it determines the most significant properties of its mechanical behaviour [66]. From the literature, the track response for this method is investigated for different scenarios based on damping property and loading type.
The steady state response due to non-oscillating moving load of a damped [67] and undamped [68][69][70] beam on elastic foundation is investigated. Other researchers studied the steady state response due to oscillating moving load. Mathews [58,59] studied the steady state response of damped and undamped beam using Fourier transformation. Duffy [71] studied the combined responses of transit and steady state of a beam on elastic foundation. Fryba [72], studied the finite and infinite models under moving oscillating load.
The second method counts for sleeper spacing and rail joints where rails are modelled as discretely supported beams. There is a variety of methods in the literature to investigate the response of discrete continuous models. Jezequel [73] analyzed infinite beam supported periodically by lateral and torsional stiffness using generalized Fourier method under non-moving oscillating load. Kisilowski et al. [74] used Jezequel's approach including a moving load on a periodically supported rail. Krzyzynski [75] modelled Eurler Bernoulli beam supported discretely using Floquet's method under a harmonic moving load. Muller et al. [76] did a comparison between the two previous methods using two beam formulations. The authors applied Generalized Fourier method using Timoshenko beam and Floquet's methods using Euler-Bernoulli beam. Smith and Wormley [77] used Fourier transform technique to do computations for one repetitive unit. The authors used periodicity condition to obtain the solutions for the other units. Resilient hinges where used by Belotserkovskiy [78] to count for rail joints. The author used Fourier transform techniques to study his model of Euler-Bernoulli beam under harmonic moving load.

Main Findings of Ballast Mechanical Behaviour using Analytical Modelling
Analytical modelling is used broadly in the literature to understand the macroscopic behaviour of railroad ballast. It is used to understand the static and dynamic responses of ballasted track under different train load cases. However, analytical modelling dose not contribute to the microscopic understanding of the mechanical behaviour of railroad ballast. Nevertheless, the use of analytical modelling enhances the knowledge, understanding and use of the defined theoretical problems that represents ballasted track behaviour.

Key Concerns and Limitations
Analytical modelling of ballasted tracks used for static and dynamic analysis. In the static analysis, ballast layer is mostly modelled as simple linear elastic spring with a constant stiffness. In the dynamic analysis, other structural elements are used to analytically model the ballast layer like masses and viscous dampers with constant damping coefficient and mass values to count for the damping and dynamic effects of ballast layer.
Analytical models of ballasted tracks have been used enormously in the literature to study the ballasted track dynamic behaviour due to their low requirement of computational time, simplicity [79] and effectiveness in the track static and dynamic analysis. Analytical modelling approach enrich the understanding of the defined theoretical problems and used as an initial, quick and simple validation to other numerical models.
Nevertheless, they have number of drawbacks in modelling the mechanical behaviour of ballast layer. While analytical modelling being useful for scoping and investigating the short-term mechanical behaviour, they lack many characteristics to analyse and visualize the complex and long-term mechanical behaviour of ballast layer.
From a material perspective, analytical approach of modelling ballast layer does not represent the real complex behaviour of railroad ballast under cyclic loading. Where, ballast behaves as elasto-plastic material with non-linear load-deformation behaviour under cyclic loading. Furthermore, the stiffness and damping properties of ballast layer is not constant with time due to ballast densifications and deteriorations under cyclic loading and after long term service. Moreover, analytical models do not have the ability to study the effect of variable parameters on ballast mechanical behaviour e.g. particle breakage, particle shape, particle size distribution and fouling.

Numerical Modelling
Numerical models of ballasted tracks are based on theoretical models that are solved numerically to describe the behaviour of each track element. With the advanced developments of computational resources, the usage of numerical modelling dramatically raises. Many researchers studied and evaluated the mechanical behaviour of ballasted track using numerical modelling approaches. There are two main numerical approaches used to model ballast layer based on material characteristics, continuum and discrete.

Origin and History
Finite Element Method (FEM) is a continuum numerical method that is widely used in many research applications. It is considered as an essential part of computer Aided Engineering (CAE). Madenci and Guven [80], Zienkiewicz and Taylor [81] and Moaveni [82] defined the FEM as a powerful computational tool that estimates the solutions of different real problems that have sophisticated domains with boundary conditions. It is introduced to solve complicated civil engineering problems related to structural and elasticity analysis. Turner [83] provided the name of FEM. The method was introduced in 1940s and its usage has increased dramatically throughout the years. FEM has been used in several real-life applications. The origin of the method was not clearly investigated. However, Hrennikoff [84] and Courant [85] were the main contributors to FEM developments. Zienkiewicz and Cheung [86] were the first people who discussed the applications of FEM by their book entitled "The Finite Element Method in Structural and Continuum Mechanics". In 1973, Strank and Fix [87] classified FEM as one of the numerical modelling methods of many physical systems. The authors accurately described the mathematical foundation of FEM which has matured and stabilized recently [88].
The first numerical models used to understand the mechanical behaviour of ballasted tracks were in the period of 1970s-1980s. Modelling the mechanical behaviour of railroad ballast has increased throughout the years as shown in Fig. 3. Although, railroad ballast is a discontinuous material, there are number of developed material models that reflect the discontinuity of ballast material. Those material models are used to model the ballast layer in FEM. Prevost and Popescu [89] introduced a constitutive model for granular material to be modelled as continuum material, despite its discontinuity properties. Wolf [90] summarized the novel approaches that used FEM to model granular material under dynamic effects. Karlsson et al. [91] used FEM to model the flow of granular material in several different geometrical silos. Kolymbas [92] discussed in his book entitled "Constitutive Modelling of Granular Materials" the different constitutive models for granular materials. Nevertheless, FEM has some drawbacks of modelling the complex behaviour of railroad ballast as discussed below.

Ballast Material Model
The main purpose of FEM is to analyze the stress-strain and deflection behaviour of continuum material under loading based on constitutive material equations. This is done by discretizing the continuum material to finite elements and solving the obtained mathematical equations of the finite elements. Each element has number of degrees of freedom and nodes based on its defined shape. The stresses and strains are numerically calculated for each node. An example of FEM of ballasted railway track consists of 33,106 frame elements and 40,517 nodes using the mesh represented in the figure is shown in Fig. 4.
The most constrained part of FEM in modelling the real complex behaviour of railroad ballast is choosing the right material model that represent the discontinuity of railroad ballast material.
Most of FEM models of ballasted tracks in the literature model the ballast layer as continuum homogenous material by discretizing the ballast layer into very small elements [93]. Other studies model ballast material as discrete connected elements using FEM [94,95]. Kuo and Huang [95] compared between two finite element models of continuum and discrete ballast layer. The first model, ballast layer was modelled as elastic continuous material using FEM. The scorned model, ballast layer was modelled as discrete connected springs using FEM. The authors found a major differences of ballast layer's settlement between the two approaches. The authors recommended that the real behaviour of ballast layer should be between the two models.
The continuum modelling of ballast layer is based on material models that reflects its behaviour. From the literature, ballast layer commonly was modelled using linear elastic material models as it requires lower computational time compared to non-linear elastic material models. Recent studies like [96][97][98] showed the importance of modelling ballast layer as a non-linear material and the authors concluded with good results of estimating the whole ballast layer behaviour. There are two methods to simulate the non-linearity of ballast. The first one which is commonly used, utilizing the Moher-Cloumb plastic principles of the ballast material model using FEM. The second one is by using hardening soil model available in PLAXIS software. Kalliainen et al. [98] and Indraratna et al. [99] used the hardening soil model in their studies about the behaviour of ballasted track. Indraratna et al. [99] recommended to use the non-linear hardening soil model in modelling ballast layer using FEM. His recommendation was based on the realistic reflection of ballast behaviour compared to large scale triaxial test  [96] results. Neverthless, Giner et al. [100] presented that generally FEM is not capable to model the short term behaviour of ballast layer under cyclic loading. Because at the initial stages of cyclic loading, the real elastoplastic behaviour of ballast is complex and difficult to be modelled using FEM where nonlinear resilient and permeant deformations occur. However, FEM is capable to model the long-term behaviour of ballast layer where the real ballast permanent deformation is minimum, and resiliency is maximum due to ballast layer densification.

Model Dimension
In the literature, the sub-structure of the ballasted track has been simulated using FEM as two, three or four layers involving the ballast, sub ballast and subgrade layers [101][102][103]. There are three types of ballasted track models using FEM based on the model dimensional; two dimensional, three dimensional and 2.5 dimensional.
In the literature, there are small number of two dimensional FEM simulations for ballasted track e.g. [104,105] as they require lower computational time compared to three and 2.5 dimensional models. The two-dimensional FEM ballasted track models do not represent the reality of ballasted tracks because of the used simplified assumptions for those models. The first simplified assumption is the plain strain condition that assumes the strain in the longitudinal direction of the ballasted track to be equal to zero. The second simplified assumption is the load application in the longitudinal direction of the ballasted track only. The previous simplified assumptions do not represent the real case in the field.
However, there are a large number of publications [98,[106][107][108] modelled the ballasted track using FEM as three dimensional. This is more reasonable as it is more representative to the real case scenario. Load distribution is considered for both longitudinal and transvers directions of the ballasted tracks Number of studies used three dimensional FEM simulations of multilayers to model the ballasted track behaviour under passing trains for number of cycles [2,94,[109][110][111]. For example, a three dimensional FEM simulations were used in modelling the ballasted tracks to understand and analyze the stress distribution of ballast layer per references [112][113][114] and ballast layer displacement per references [98,107,115] from a macroscopic point view.
Other researchers modelled the railway track as 2.5D to reduce simulation computational time [116][117][118][119][120]. The 2.5D is used to model the geometry of the track as two dimensional while considering for three-dimensional load cases. Figure 5 shows examples of different dimensional ballasted track models using FEM.

Main Findings of Ballast Mechanical Behaviour using FEM
FEM is a powerful tool in modelling the entire ballasted track behaviour and analyzing the track behaviour for large number of loading cycles. Moreover, FEM is an efficient numerical tool in the dynamic analysis of entire ballasted track. It can study and analyze the macroscopic behaviour of ballast layer with response to static and dynamic loading. However, the microscopic behaviour of ballast can't be analysed and visualized using FEM. For example, force chain distribution, particle displacement and particle breakage. Additionally, some parameters that influence the strength of ballast layer can't be studied using FEM because ballast layer is modelled as a continuum material. For instance, particle shape, particle size distribution and fouling. From the literature, most of FEM models of ballasted tracks were static with using dynamic amplification factors to consider the dynamic affects. There are three types of FEM models of railway track that consider the dynamic effects [100]: complete FEM model of railway track like [127,128]; FEM model of railway track where finite elements applied in the model's boundaries like [129,130]; FEM model of railway track where boundary elements applied in the model's boundaries like [131][132][133][134].
Most of FEM models of ballasted tracks used same stiffness parameters for static and dynamic cases in their models [100]. Ganesh and Sujatha [115] analysed the ballasted track modulus using FEM for different types of sleepers (pre-stressed concrete and wooden sleepers) and soil conditions. The authors showed the significant contribution of rail pad and ballast stiffness parameters to track displacement.
Other publications in the literature studied the influence of train speed on the mechanical behaviour of ballasted tracks using FEM. El Kacimi et al. [135] found that at track critical velocity, ballast track deformation is maximum. Varandas et al. [114] found that at track critical velocity, the stress path is significantly influenced, and principle stress rotation of ballast is maximum.

Key Concerns and Limitations
FEM is useful tool in modelling the long-term overall ballasted track behaviour under different loading conditions and large number of cycles. It provides the macroscopic behaviour of railroad ballast layer. Few FEM models used to validate experimental test results like box test [48,[136][137][138] and direct shear test [139]. The requirement of computational time is lower compared to Discrete Element method.
The most significant limitation of FEM in modelling the mechanical behaviour of railroad ballast is the application of the material model that reflects the realistic discontinuous and elastoplastic properties of railroad ballast. Although the continuum assumption works in most cases, the insight visualization of stresses and displacement cannot be correctly evaluated using FEM because ballast layer is modelled as a continuum material. For instance, it is difficult to model ballast layer Initial settlement, volumetric change, particle breakage, force chain distribution, particle displacement, particle shape, particle size distribution and fouling using FEM.

Origin and History
Discrete Element Method (DEM) is a numerical method used to solve mathematical problems associated to discrete characteristic material like granular material. Each particle has its own properties of displacement, velocity, acceleration and contact forces. The calculation process for each property of each particle in a granular assembly is a tremendous complex. Therefore, discrete characteristics of granular material make the understanding of its mechanical behaviour very difficult. DEM is a powerful tool that can analyze the mechanical behaviour of granule material both microscopically and macroscopically [140]. There are two approaches used in DEM based on particle contact nature; hard and soft spheres.

Hard Sphere DEM
The hard sphere method considered the contact between particles to be rigid (Fig. 6). The overlaps and deformations of particles are not simulated using this approach. Forces between particles are impulsive and implicitly considered. Only momentum exchange is happened due to collisions. Particles' velocities are explicitly calculated based on the material properties (coefficient of restitutions). One collision at time is considered. Never two or three collisions are happened simultaneously. Typical application of hard sphere method is rapid granular flow simulations [141]. The origin of hard sphere DEM is not clearly identified. It is stated that Metropolis et al. [142] introduced the hard particle system models in the study of thermodynamics equilibrium of hard disks system using a Monte Carlo method. From their publication, Molecular Dynamics (MD) theory developed. Alder and Wainwright [143] introduced the theory of DEM by their molecular dynamics studies. As well as the shear viscosity studies that were done on a liquid argon using MD by Ashurst and Hoover [144]. MD calculates the forces based on the distance between particles and it considers the particles to be frictionless with no rotation.
In the early 1980s, application of MD for large particles is introduced as a Granular dynamics (GD) to model rapid granular flows by [145][146][147][148]. GD is considered as the outgrowth of hard sphere DEM type. It considers the particles to be rigid. The collisions are not simultaneously simulated. It is known as event driven simulation. At each time step one collision (contact) is considered. Continues or simultaneous contacts are not considered. This method is applicable for granular systems that has quick contact interaction between particles like rapid granular flows [149]. Forces are calculated implicitly based on particle material like surface friction and coefficients of restitution.

Soft Sphere DEM
The soft sphere method considered the contact between particles to be rigid partially. The soft sphere approach considers the overlaps and deformations of particles during the simulation [130]. Soft sphere method can handle number of particle contacts [129]. Soft sphere DEM has been used broadly to study many discrete particles phenomena. From the literature, most of the DEM models for railroad ballast used soft sphere discrete element method. Few of them used hard sphere discrete element method like [131,132]. O'Sullivan [133] drew a schematic diagram to illustrate and summarize the differences between hard and soft sphere DEM approaches as shown in Fig. 7. Cundall and Hart [150] set the basics of soft sphere DEM. The authors defined DEM as a calculation tool that recognizes the contacts and allows the rotations and displacements of each discrete element. In 1971, the first computer model that models the progressive failure of a rocky discrete elements was introduced by [151]. The model was based on the friction and normal stiffness for the interaction between discrete blocks.
In 1979, the distinct element method which is commonly known now as DEM was introduced to model granular assemblies including particles overlapping by [140]. The purpose was to simulate the particle contact force chain in the assembly. In his computer program "Ball", he utilized the aspects of calculation cycle, law of force displacement, law of motion and damping effects to simulate the contact force distribution. The author studied the internal mechanisms of 100 and 1000 discs and their responses to stress. The initial state of 100 discs is shown in Fig. 7. Initial contact force chain in 100 discs by [140] The number of publications related to modelling railroad ballast using DEM has been increased throughout the years as shown in Fig. 8. It became a required tool for many engineering applications, especially in 1990s when computers reached high level of performance. Soft sphere DEM has been used broadly in many research areas and DEM software packages. In railway ballast applications, most of the DEM models for railroad ballast in the literature used the soft sphere discrete element method. Saussine et al. [152] and Lim and McDowell [153] were the first people who simulate the mechanical behaviour of railroad ballast using DEM.

Ballast Contact Model
In DEM there are two types of interactions between the distinct elements; particle to particle and particle to wall. The main purpose of contact models is to simulate both interactions during the simulation period based on theoretical knowledge of contact mechanics science. Modelling the contact mechanics of granular materials using DEM has been discussed extensively by geotechnical scientists [102][103][104]. The contact models used in DEM of railroad ballast are discussed in this section.
The contact between two particles is not at a single point but on a very small finite area because of the particle's deformation. At the contact area there are normal and tangential contact plans. Therefore, total contact force consists of normal and tangential forces. It is not easy to accurately simulate the contact area between particles due to many factors like geometrical shape and material characteristics [141].
The accuracy of DEM results is dependable on the accuracy of material properties parameters used in the contact model [141,154]. DEM commonly uses simplified contact models that find the forces and torques acting on a particle due to contact, in the sake of efficient computational time. There are different contact models for different element shapes (i.e. spheres, polyhedrons or others). Most of the contact models are developed for spherical contacts using springs and dashpots [141,155]. Mishra [156] discussed the different contact models used for spherical shapes. However, the most used contact models for railroad ballast modelled as spherical shapes are linear elastic contact model by [140] and non-linear elastic Hertz-Mindlin contact model by [157,158].
The linear elastic model is simpler and requires less computational time compared to the non-linear elastic Hertz-Mindlin contact model. There were many comparative studies e.g. [159,160] between the linear elastic and non-linear elastic Hertz-Mindlin contact models and the authors found that results using both models are close and in good agreement. However, Di Renzo and Maio [159] recommended the Hertz-Mindlin contact model to be used in modelling granular behaviour using spherical shapes that can provide an accurate and deep understanding to granular motion.
Non-spherical shapes like polyhedrons have their own complicated and time-consuming contact models that are based on mathematical equations. For polyhedrons, there are few methods to calculate and detect the contacts between them. For example, common plane method [161,162], overlapped volume intersection by [163] and potential particle shapes by [164]. Table 1, summarizes the contact models used in the literature to model railroad ballast. Obtained from the Scopus using the following keywords: railroad ballast OR ballast track OR ballasted track AND discrete element model OR discrete element method OR DEM

Particle Shapes
There are different particle shapes used to represent railroad ballast in DEM. For realistic results, the modelled distinct element should have almost the same real shape. The complexity of particle shape in DEM is a trade-off with computational time and result's accuracy. The more the particle shape used in DEM model is complex, the more the computational time is needed to detect the contacts, but DEM results are more realistic.
Particle shape should be efficient to accurately describe the real shape and at the same it should be simple enough to reduce the computational time. Many studies e.g. [165][166][167][168] presented the importance of particle shape of railroad ballast used in DEM simulations.
Spheres in 3D or circles in 2D are the mostly used shapes in DEM applications. As they require low computational time and most of contact modelled are developed for spherical shapes as discussed above. Lim and McDowell [153] and Lobo-Guerrero and Vallejo [169] introduced the DEM modelling of railroad ballast using spheres and disks respectively. However, spheres and disks do not reflect the real behaviour of railroad ballast due to their easy rotation under loading and they do not provide the interlocking property of railroad ballast [153,166,169,170]. Therefore, three approaches have done in the literature to model railroad ballast with taking into consideration the shape effects.
The first approach is multi-sphere. Its main advantage is that it almost considers the irregularity and angularity of ballast meanwhile it sustains the computational efficiency. It gives a better approximation to the real ballast shape than spheres. It is used in many other studies like geotechnical [171][172][173] and aerospace [170,174]. In this approach, the railroad ballast shape is approximated by a number of overlapping or touching spheres to form a cluster or sometimes identified as clumps (Fig. 9). The bond between the spheres is normally set to infinity to form a rigid cluster. Particle breakage in this case is not allowed. The contact forces calculation is done only between clusters not spheres in each cluster. Some researchers set a finite bond strength between the spheres to model particle breakage and crushing [169,175,176].
Lu and McDowell [166] developed simple steps to generate clumps. Their technique aims to overlap different ball numbers and diameters to form complex clumps. It includes the aspects of ballast angularity, sphericity, and surface roughness. Mahmoud et al. [5] developed novel two multicircles approaches that capture the angularity of railroad ballast using MATLAB and AutoCAD routines in addition to the typical routine of spherical clumps. MATLAB routine generates hexagonal close-packed assemblies. A closed hexagonal shape is used for each scanned ballast shape and then it is filled by the same sized circular elements. AutoCAD routine represented by filling each scanned ballast with number of tangential different size circles. The two innovative routines allow the modelling of particle breakage. The 2D scanned ballast images were taken from Aggregate Imaging System (AIMS) of the authors' database. Another new variation of multi-sphere approach is introduced by [177][178][179]. The authors used clumps with bonded asperities to model the particle abrasion as shown in Fig. 9. Chen et al. [53] used simplified clump shape of two overlapped spheres to model the large number of particles and load cycles with low computational time. Zhang et al. [180] introduced a new approach to generate multi spheres clusters of bonded spheres using laser scanning technology. First, point cloud data of the ballast surface is taken by laser scanner. Second, based on the point cloud a closed surface is made to model particle geometry. Then, spheres are generated within the cloud surface. Figure 10, summarizes the flow diagram of their approach. Indraratna et al. [181] recommended the proper number of spheres in a cluster to be 10-20 spheres. Many studies [153,166,168,178] showed the realistic deformation behaviour of railroad ballast modelled as clusters compared to spheres this is due to interlocking property of clusters.  [20,198,204,211,212,214,215,223,224] Other contact models based on mathematical equations (non-spherical shape) [33,37,41,163,164,183,185,190,196,206,[225][226][227] Fig. 9 Model of ballast particle as a sphere, b cubical cluster of 8 spheres, c tetrahedral clump of 10 un-breakable spheres with 8 breakable asperities [178] The second approach is to model the ballast shape as polyhedrons. A polyhedron shape is defined by the number of corners, edges and faces (Fig. 11). Its main advantage is that it represents the realistic shape of ballast. However, it is very time expensive due to the massive time needed to detect and calculate the contacts of each particle, especially edge contact force [170]. Alonso-Marroquín and Wang [182] presented a new algorithm to speed up the computational time by deleting the edges of polyhedrons. Another disadvantage of using polyhedrons shapes in DEM models is that there are low numbers of descriptive contact models for this shape. Contact detection algorithm and contact force calculations should be designed. Most DEM software use sphere shapes. It can be used for low number of particles and load.
Saussine et al. [183] modelled the ballast behaviour under cyclic loading. The main objective of their study was to show and validate the ability of DEM to model railway ballast. Ballast particles were modeled as pentagonal shape. Eliáš [163] studied the behaviour of polyhedron modelled railroad ballast under monotonic loading. Small number of particles are modelled due to computational time issues. The author used a novel approach to generate random polyhedral particles using Vorni tessellation technique with three geometrical ratios as shown in Fig. 12. The results are promising compared to experimental work. Ahmed et al. [164] simulated the triaxial test of railroad ballast using polyhedron shape.
The third approach is to introduce a rolling resistance moment at the particle's contacts. This approach was introduced by Jiang et al. [184]. The main idea of this approach is to add a rolling resistance moment in the spherical particle contact to resist the rolling motion. This represents the angularity and interlock properties of real ballast. This approach is simpler and requires lower computational time compared to the previous two approaches [185]. The contact detection and contact force calculation are straight forward as the particle consists of one sphere with known radius centre and position. The number of spheres is lower with relative to multi-sphere approach. Therefore, less computational time is needed to detect the contacts and calculate the contact forces. There are many studies not related to railway engineering which used this approach [186,187]. Irazábal et al. [185] recently used this approach with an application to

Fig. 11
Real ballast and corresponding polyhedron simulated ballast [164] railway engineering to study the lateral resistance of railroad ballast using DEM. Table 2, compares between the used particle shapes in modelling railroad ballast using DEM. Table 3, summarizes the used particle shapes in modelling railroad ballast using DEM in the literature.

Simulated Train Loading Type and Loading Cycle Number
The train consists of a number of cars. Each car has typically four axles with different spacing. Each axle exerts a load on ballast layer. The loading from the train depends on different    [6, 20, 38, 53, 165, 166, 175-181, 191, 198, 210, 214, 217, 220-223, 228, 229] Polyhedrons or polygons [33,37,41,163,164,169,180,190,194,195,204,219,226,227] parameters. For instance, car length, car weight, axle spacing and time between passing trains. From the literature, most of the studies related to railroad ballast modelling using DEM simulate train loading as pure continuous loading as shown in Table 4. There are a few studies which simulated the behaviour of railroad ballast under haversine loading. However, as recommended by [188] that haversine can represent one single axle loading only. A number of loading axles cannot be represented by haversine [7].
From the literature, DEM was used to understand the short-term behaviour of railroad ballast under small number of loading cycles. Long term behaviour of ballast behaviour after a large number of loading cycles is not yet understood using DEM. The used number of loading cycles using DEM is very low compared to experimental work. Ngo et al. [6] and Chen et al. [53] used 500,000 and 200,000 loading cycles in their experimental tests respectively. However, the maximum number of loading cycles used in DEM to understand the mechanical behaviour of railroad ballast is 6000 for a two dimensional test model and 4000 for a three dimensional test model by [175] and [6] respectively, because DEM requires huge computational time.
Computational time of a DEM model limits the researchers from various perspectives to be considered in the modelling process. For example, the number of loading cycles, particle shape, particle breakage, simulation dimension and loading type. It is very difficult to simulate ballast behaviour using DEM under a large number of loading cycles using three-dimensional scale and complex particle shapes including particle breakage. Some studies used large number of loading cycles but in two-dimensional scale. Others used complex particle shapes and included particle breakage but for a low number of loading cycles. It is difficult to consider all the perspectives of realistic modelling of railroad ballast in one DEM simulation as shown in Table 4. It is a trade-off between realistic simulation and computational time.

Software Packages
Cundall and Strack [140] are the first people who computationally implemented DEM with a code named "BALL". Then, a variety of DEM packages have been used to computationally implement DEM with advanced criteria including huge number of distinct particles in different loading and environmental conditions. Table 5 below, describes the used DEM software packages in the literature to model the mechanical behaviour of railroad ballast.

Shear Strength
Ballast shear strength is an important criterion in selecting ballast material. Ballast shear strength is gained from the particle angularity and texture of railroad ballast. Ballast shear strength preserve the functionality of railroad ballast in long term behaviour. It resists the slip failure of ballast which could cause track realignment and train derailment [4]. Many studies have studied ballast shear strength via numerical DEM modelling to validate experimental tests or to study the effect of different mechanical, physical and geometrical properties on the  [189]. The authors modelled the large-scale triaxial test using PFC and analyzed the volumetric change under triaxial conditions. The authors observed a dilation to the modelled ballast layer's volume rather than contraction at high confining pressure which contradict with the experimental results. The authors explained this due to not considering particle breakage in their model. However, Lu and McDowell [178] showed more realistic volumetric change behaviour by modelling railroad ballast as clump with breakable bonded asperities. The authors studied the effect of different parameters on the ballast mechanical behaviour under monotonic triaxial loading. The parameters are particle shape, angularity, interlocking, abrasion and friction coefficient under a range of confining pressures. The authors analyzed both the deviator and volumetric changes under a slow shearing rate by the application of two walls (top and bottom). Which was later emphasized by [177]. Qian et al. [190] studied the effects of slow and fast shearing rate on ballast behaviour using DEM. The authors found that shearing rate has a negligible effect on the ballast shear strength. Indraratna et al. [191] studied the influence of load frequency and particle breakage on the ballast mechanical behaviour under cyclic loading through DEM. The authors found that load frequency influenced the permanent deformation of ballast. Ngo et al. [38] used DEM to model the microscopic behaviour of clay fouled ballast. The authors observed using DEM that number of contacts increased significantly as level of fouling increased but the maximum contact force decreased due to clay fouling. The authors emphasized on the ability of modelling granular material using DEM which allow the using of different particle sizes. The authors concluded the significance of modelling particle size as they influence the bulk shear strength of granular material. Huang et al. [192] were the first people who used DEM to model direct shear test of fouled and clean railroad ballast based on their experimental work [17]. The authors analyzed the stress-strain behaviour due to direct shearing of fouled ballast using DEM. However, the volumetric change of fouled and clean ballast using DEM model was not investigated. Indraratna et al. [181] studied the stress-strain, volumetric change and contact force distribution of clean and fouled railroad ballast under direct shearing using DEM. The authors found that the stress-strain behaviour of numerical DEM model in good agreement with experimental results. However, volumetric change behaviour of numerical model showed a higher dilation behaviour under direct shearing compared to experimental results. This was due to different particle shape angularity between the modelled particle shape and real ballast shape.

Direct Shear Test Using DEM
Wang et al. [20] analyzed the stress-strain behaviour of ballast during direct shearing under different normal stresses (100,200 and 300 kPa). The authors have concluded that the shear strength increases with the application of normal stress. Their work included the effects of particle size distribution, particle shape and normal stress on the shear strength of ballast using DEM model of direct shear test. The authors recommended particle size distribution in the model should be the same to the particle size distribution in the experiment. The authors concluded that railroad ballast should have broad grading that increase the shear strength. Their conclusion was the same of Indraratna et al. [8,193]. The authors did not analyze the volumetric variation in their study as Indraratna et al. [180] did.  [5, 6, 20, 38, 53, 153, 165, 166, 169, 175-178, 181, 191, 198, 210, 211, 216, 217, 220-222, 229, 231] EDEM by Dem Solutions [215,223,224] In house BLOCKS3D by the University of Illinois [232][233][234] [ 33,37,41,180,190,195,219,226,227] Interactive graphical software to create ballast particles developed by Ahmed et al. [164] [ 164] Zaho et al. [194] showed the importance of particle angularity in modelling of ballast behaviour. The authors studied the shear stress-strain, volumetric change and the relationship of friction angle and cohesion to angularity. The main conclusion of their study was that polyhedral particles gives more realistic behaviour of real ballast under shearing unlike spherical particles. However, the authors mentioned that the massive computational time was a drawback for their study.
Other studies modeled direct shear test to calibrate material properties of fouled and clean ballast with relative to experimental work [195,196]. The authors validated their model using Bulk Calibration method. By trial and error, the authors have changed some material parameters like friction angle until the numerical direct shear DEM model results match the direct shear experimental results. Huang and Tutumluer [195] used the validated parameters from direct shear DEM model to study the displacement behaviour of fouled and clean ballast using a half-track DEM model. Moreover, Huang and Tutumluer [196] used direct shear test to validate their novel approach of particle generation in DEM that considers the particle shape effects. Their novel method called image-aided DEM that used BLOKS3D software.
Some studies used DEM numerical models of direct and triaxial shear tests to show the advantages of using geogrid reinforcement in enhancing ballast shear strength. Qian et al. [33] presented a comparative study between the different shapes types of geogrids used to reinforce railroad ballast. Miao et al. [165]. Visualized comprehensively the microscopic behaviour of ballast-geogrid interaction by modelling the pull-out behaviour of a triaxial geogrid inserted in a ballast layer using DEM. Chen et al. [53] showed the ability of geogrid on reducing lateral and vertical displacement under different cyclic and confining loadings through Process Simulation Test (PST).

Ballast Breakage
After long term service, ballast particles break. There are two types of ballast breakage. Corner breakage and particle splitting breakage. Lackenby et al. [197] found that corner breakage occurred mostly under dilation condition (low confining pressure) and splitting across particle breakage occurred mostly under contraction condition (high confining pressure).
Many researches presented the importance of considering ballast breakage in DEM models [169,191]. However, several researches related to railroad ballast did not include particle breakage in their DEM model due to the high requirement of computational time [153,166]. Table 6, summarizes the literature about the use of particle breakage model in modelling the mechanical behaviour railroad ballast using DEM.
Ballast breakage is modelled in DEM using two approaches. The first approach is to identify a finite bond tensile strength between the particles of the clumps. The breakage failure occurs on the bond between the particles in each clump; when the acting contact stress is more than the bond tensile strength, bond breakage occurs. This approach is used when ballast is modelled as spherical clumps. This approach was presented by [176,180,191]. Zhang et al. [180] used a complex particle shape with bonded spheres to represent railroad ballast shape and breakage respectively. This approach is very time consuming. The authors studied the effect of bond strength (10 and 30 MPa) used in DEM model on the permanent deformation of ballast. The authors visualized that most of the broken bonds were under the sleeper. The authors found that low bond strength between spheres in one clump produced additional settlement to ballast layer. Some studies [179,198] used box test to investigate ballast abrasion (corner breakage). Ballast particles were simulated as unbreakable clumps consists of large particles with small bonded spheres (asperities) attached to the large particles. However, this approach failed to model particle splitting breakage. But it succeeds in modelling the ballast abrasion (corner breakage).
The second approach is to identify a tensile strength (critical stress) of the clump by which after this tensile stress the cluster original geometry is replaced by different smaller sized fragments. This was presented by different researchers using circular clumps [5,169,175,199] and polyhedron [163] as shown in Figs. 13 and 14 respectively.

Effect of Confinement
Several studies were conducted to understand the effect confining pressures on the mechanical behavior of railroad ballast using DEM.
Some researchers studied the effect of confinement on ballast particle breakage experimentally [24,25] and numerically using DEM [175]. The conclusions of the previous studies agreed that confinement influenced ballast particle breakage. High confinement reduces ballast particle breakage. At low confining pressure (below 10 kPa), ballast breakage is maximum.
Other studies showed the significant of confinement on the lateral resistance of ballast layer. It is recommended to provide ballast layer with good confinement pressure to reduce ballast lateral and vertical displacements [24,25,53,200,201].
Although confining pressure is critical in designing railway tracks, there is no standards exist that relates to the  [175,179,180,198] Not including particle breakage [6,53,153,166,204] confining pressure of railway tracks [25]. Selig and Alva-Hurtado [202] reported that the provided confining pressure to the ballast in the field is in the range of 5-40 kPa, which is small relative to applied vertical loadings from the trains. Indraratna et al. [24] recommended to use a confining pressure in the range of 30-90 kPa for a cyclic deviatoric stress of 500 kPa. Lackenby et al. [25] recommended almost the same range of confining pressure 25-95 kPa for a cyclic deviator stress of 500 kPa.

Settlement Behaviour
The most benefits of numerical modelling using DEM that particle displacement and force chain distribution can be visualized unlike the Finite Element method. This section describes how researchers used DEM to understand the ballast displacement through modelling box test, half-track model and full track model. Modelling the full track or half-track using one or more sleepers is very computationally expensive. Therefore, modelling box test were commonly used in the literature to study the effect of geometrical and physical ballast properties on ballast mechanical behaviour using DEM meanwhile preserving the computational efficiency.
The conclusions of the studies regarding ballast settlement behaviour are almost similar to experimental results that ballast particles beneath the sleeper displace downward and beside the sleeper displace upwards. However, some of the studies include in their simulations different parameters related to railroad ballast to investigate their influence on ballast settlement as summarized in Table 7. For instance, ballast particle shape, ballast fouling, geotextile reinforcement, under sleeper pads, loading frequency and confining pressure.
Lim and MacDowell [153] and Lu and MacDowell [179] were the first researchers who used DEM to understand the mechanical behaviour of railroad ballast through box test. Their main concern was on validating their DEM model and Fig. 13 The induced tensile stress on a simulated aggregate and the resulted shape after breakage [169] Fig. 14 The induced particle stresses and resulted aggregate shape by [163] [198] Influence of breakage on ballast vertical settlement [169] Influence of train speed on ballast vertical settlement [205] Influence of particle size distribution on ballast vertical settlement [206] showing the importance of DEM in visualizing the microscopic behaviour of railroad ballast. Lu and MacDowell [166] discussed the effect of geometrical shapes on ballast settlement and rotation under sinusoidal cyclic loading. The authors were the first people who showed how DEM can be utilized to visualize the railroad ballast particle displacement and rotation. Ballast displacement was shown as vectors and the thickness of the vector is proportional to the displacement magnitude of the particle. The authors pointed out that ballast displacement behaviour using clumps is almost similar to the box test results obtained by [203]. Hossain et al. [175] discussed ballast settlement using DEM through box test model under sinusoidal cyclic loading and different confining pressures. The authors emphasized on the significance of confining pressure in reducing ballast settlement and concluded that around 2000 cycles ballast settlement is maximum at all confining pressures. The same conclusion was observed by [53].
Ngo et al. [6] studied the influence of geogrid on fresh and fouled ballast settlement behaviour; experimentally using process simulation test and numerically using DEM under cyclic loading. The authors observed that geogrid reduced ballast vertical and lateral settlement for all fouling indices. Li and McDowell [198] studied the influence of under sleeper mat on ballast settlement behaviour using DEM through box test. The authors found that under sleeper mat reduced ballast settlement and particle abrasion under sleeper.
Zhang et al. [180] presented the settlement behaviour under simulated train loads of heavy haul vehicle passage through DEM model of box test. Ballast was modelled using complex clumps with bonded spheres to represent ballast shape and breakage respectively. This approach is very time consuming as they include large number of particles in one simulation. Furthermore, including particle breakage is an additional calculation to the DEM simulator. The authors showed the effect of modelling particle breakage using different bond strengths (10 and 30 MPa) on the ballast settlement behaviour. The maximum sleeper settlement for one haul car passage was in the range of 0.83-1 mm.
Ji et al. [204] investigated the influence of cyclic loading frequency on ballast settlement behaviour. The authors found that cyclic load frequency affects ballast settlement behaviour. Maximum settlement was observed for high load frequency.
Other researchers prefer to model the ballast settlement behaviour under one sleeper or number of sleepers to reflect the track/half-track condition. Lobo-Guerrero and Vallejo [169] discussed the influence of modelling particle breakage on settlement behaviour of ballast. The authors did two 2D DEM models of track section consists of three sleepers as shown in Fig. 15. The loadings are applied at the same time on the three sleepers which is not the real case at the field. The authors found that although few particles were crushed in the model, ballast displacement increased extremely.
Huang and Chrismer [205] studied the effect of train speed on ballast settlement. The authors developed a 3D model consists of 3 sleepers at the top of ballast as shown in Fig. 16. The authors studied two cyclic simulated loading cases: The first one for freight train and the second one for mixed traffic (freight and high-speed trains). The authors found that train speed influenced the track responses. High Fig. 15 The simulated track model in PFC 2D by [169] Fig. 16 Three-dimensional half-track DEM model [205] speed trains enlarged track responses and introduce some uplift track displacement in front of the train.
Bian et al. [206] discussed the effect of particle size distribution using different ballast standards on ballast settlement using half-track DEM model (Fig. 17). The authors have used the standard ballast gradations for three countries Australia (Rail Infrastructure Corporation and Queensland), France and United States (AREMA No. 24, No. 3, and No. 4 gradations). The applied cyclic load was from dynamic track model done by [207]. The authors concluded that gradation No.24 by AREAM resulted the least track settlement under loading relative to others.
Zhang et al. analyzed the displacement behaviour of ballast layer under 10 sleepers as shown in Fig. 18. The authors concluded that the effected ballast zone is within seven sleepers under the bogie. Ballast under the loads directly displaced downwards, but ballast beside the applied loads displace longitudinally and vertically.

Key Concerns and Limitations
DEM is a powerful tool for understanding and visualizing the mechanical behaviour of railroad ballast. DEM can visualize the microscopic and macroscopic behaviour of railroad ballast through a comprehensive insight into the contact force distribution and particle displacement during loading. DEM can study the effect of different parameters of railroad ballast like particle shape, particle size distribution, particle breakage and fouling on the mechanical behaviour of railroad ballast.
Many Researchers used DEM to understand the mechanical behaviour of railway ballast. Most of the DEM models for railroad ballast were used to calibrate material properties or validate experimental tests. For example, particle crushing test [153,208,209], uniaxial test [153,163,210,211], triaxial test [37,41,177,178,191,208,212], direct shear test [20, 181, 192, 194-196, 213, 214] and box test [6,53,153,166,175,179,180,198,204,215,216]. Some researchers modelled ballast layer under number of sleepers [5,169,205,217,218] and under one sleeper [183,195,206,219] to represent the track/half-track condition using DEM for a low number of loading cycles.
However, Full track analysis under dynamic train loading for a large number of loading cycles is not fully understood using DEM. The main limitation of DEM is the massive requirement of computational time despite the advancement of computational resources. The effect of dynamic loading due to track and rail irregularities on the mechanical behaviour of railroad ballast has not been fully investigated yet. Most of the studies on modelling the mechanical behaviour of railroad ballast modelled ballast under quasit-static loading.
The length scale and dimension of the DEM model and number of loading cycles are a trade-off between realistic modelling and computational time. Large-scale and threedimensional ballasted track DEM model using large number of loading cycles reflect the real ballasted track condition. Nevertheless, it is very difficult to model this case using DEM as it requires huge computational time.
From the literature, DEM was used mainly to validate experimental tests related to railroad ballast for a small-scale   [217] model. Most of the studies modelled the behaviour of railroad ballast for short term using a low number of loading cycles under pure continuous sinusoidal loadings. Table 8 below, summarizes the differences between analytical methods, FEM and DEM in modelling railroad ballast behaviour.

Conclusion
The conventional ballasted tracks are the first introduced type of railway tracks to the world. Ballasted tracks have been used widely worldwide. Ballasted track consists of two main structure; superstructure and substructure. Superstructure consists of rails, rail pads, sleepers and fastening system. Substructure consists of ballast, sub ballast and subgrade. Ballast is considered as the basic element in ballasted track due to its significant contributions to track stability. Ballast maintenance is needed periodically after service to ensure the functionality of ballast layer. Ballast maintenance needs huge financial support beside that it is a time-consuming activity. This raise the research work on understanding the mechanical behaviour of railroad ballast to design and maintain ballast layer effectively. From the literature, numerical approach has been used intensively to understand and estimate the mechanical behaviour of railroad ballast. With the developments of computer resources numerical methods like FEM and DEM have been used broadly to simulate ballast behavior.
FEM can model the entire track behaviour with an efficient computational time. FEM is a good numerical tool to explore the macroscopic behaviour if ballasted track. The challenge of FEM is to choose an appropriate material model that simulate the complex discontinues ballast material behaviour.
DEM is alternatively used, where discontinues property is taken into consideration. DEM can visualize the microscopic mechanical behaviour of ballast under loading. The initial settlement, volumetric change, particle breakage, force chain distribution, particle displacement, particle shape, particle size distribution and ballast fouling can be investigated using DEM.
From the literature, DEM was used mainly to validate experimental tests related to railroad ballast for a smallscale model. Most of the studies modelled the behaviour of railroad ballast for short term using low number of load cycles under pure continuous sinusoidal loadings.
It is necessary to consider the following for future research work related to geomodelling of ballast behaviour; full track DEM analysis under a more realistic dynamic train loading for a large number of loading cycles.