A discrete element study on the deformation and degradation of coal-fouled ballast

This paper presents the results of Discrete Element Modelling (DEM) which quantitively examine the effect of coal fouling on the deformation and degradation of ballast upon cyclic loading. The degradation model described herein considers the Weibull distribution effects in tandem with a granular medium hardening law that incorporates the maximum contact criterion to capture surface abrasion and corner breakage of angular ballast. The DEM model had been calibrated initially with laboratory data obtained from large-scale direct shear testing. Subsequently, a series of cubical shear test simulations have been carried out using DEM to understand the behaviour of fouled ballast whereby the numerical particle degradation modelling could simulate the experimental response of the ballast assembly at various fouling levels. The results show that the increased level of fouling exacerbates the sleeper settlement, while decreasing the resilient modulus and the particle breakage. Ballast beneath the sleeper experiences significant breakage compared to the crib ballast, and not surprisingly, the extent of damage decreases with depth. Rigorous microscopic analysis is also presented in relation to inter-particle contacts, particle velocity and anisotropy of the ballast assembly. This micromechanical examination highlights that the decrease in ballast breakage for fouled assemblies is predominantly attributed to the inevitable decrease in inter-particle contact pressures as effected by the coating of ballast aggregates by the coal fines.


Introduction
Ballast consists of coarse-grained particles that are commonly placed under sleepers as a load-bearing stratum. The main functions of ballast are (1) distributing the high imposed sleeper/ballast pressures to subballast (capping) and subgrade layers at a reduced stress level; (2) maintaining vertical and lateral alignment of tracks; and (3) providing sufficient permeability for free drainage of track [18]. During train loading, ballast aggregates undergo continuous degradation and breakage and as a result, ballast becomes fouled [45]. The fouling is also induced by the infiltration of fines, for instance, pumping of fluidized subgrade and coal dropping from carriages to the ballast track [12,42,50]. The ballast fouling has been reported as the primary source to many hazardous track conditions, such as impeded drainage, differential settlements at track sections, and reduced load-bearing capacity of ballast layer [13,48].
Ballast breakage caused by cyclic loadings has been studied by means of laboratory tests, and results have shown that the breakage is governed by characteristics of individual particles (e.g. shapes, sizes and angularity of particles), loading and boundary conditions [11,20,22,49]. However, most of the experimental research has been carried out in a way that representative model elements where in-situ aggregate is prepared and subjected to simplified stress conditions (e.g. confining pressures in triaxial tests, normal pressures in direct shear tests) along predetermined failure planes. Research studies have also been carried out to capture the dynamic responses of ballast aggregates under a realistic rail track condition [28,32]. Large-scale cubical box tests which accommodate a section of sleeper and ballast with sufficient depths ( ! 300 mm) have been carried out to study the load-deformation responses of ballast in recent times [14]. The mechanical performance of ballast subjected to monotonic and cyclic train loading can be well studied with ballast/sleeper interaction being considered more intuitively. Despite significant progresses that have been obtained so far, most of the existing research on ballast deformation and degradation is limited at the macroscopic scale due to the limitation of laboratory tests in capturing the microscopic response of ballast.
The discrete element method (DEM) pioneered by Cundall and Strack [8] provides a powerful tool to explore the performance of granular materials from a particulate perspective [4, 10,25,29,43]. By adopting the DEM, Lim and McDowell [26] and Lobo-Guerrero and Vallejo [31] numerically studied the effect of ballast degradation on permanent deformation in the cubical box test with satisfactory results, albeit the simulated ballast particles are in the form of assumed simplified shapes. Indraratna et al. [16] and Ngo and Indraratna [42] also applied the DEM to investigate the degradation behaviour of ballast in the presence of reinforced geogrids or under-sleeper rubber pads. In the DEM simulation, the breakage of ballast is normally simulated by either replacing the broken particles with several self-similar fragments while maintaining the same volume [31] or removing the attached spheres (3D) or cylinders (2D) from their parent clumps after reaching the maximum bond strength [26,42,52]. However, the existing particle degradation models cannot well reflect the continuous breakage process that happened to ballast particles in the real track fields. Therefore, a particle degradation model that is capable of capturing the continued abrasion and breakage of ballast particles is required for a better understanding of the dynamic responses of ballast aggregates, as well as the underlying mechanisms of particle breakage.
In the presence of fouling materials, the fabric characteristics of ballast are continuously changed (i.e. due to ballast breakage and accumulation of external fines), which had been investigated to adversely impede hydraulic conductivity [51] and to reduce bearing capacity [50,9] of the ballast layer. Laboratory studies carried out by Indraratna et al. [14] and Tennakoon and Indraratna [50] have demonstrated that the clogged fines also play a significant role in ballast degradation. The detrimental effect of fouling materials on ballast has also been investigated by numerical simulations [12,43,46], albeit the degradation of ballast particles has not been considered in those studies. Ngo and Indraratna [41] have performed DEM simulations to explore the degradation behaviour of fouled ballast. However, the simulations were carried out on triaxial specimens under monotonic loading conditions. The cubical shear box tests under cyclic loadings were simulated by Ngo and Indraratna [42], but the effect of fouling materials was not considered in their study. Indraratna et al. [14] modelled the cubical shear box tests on fouled ballast under cyclic loading conditions adopting a 2D particle degradation model, but the corner abrasion and surface attrition of ballast particles were not well captured.
Therefore, in this study, a 3D particle degradation model that can capture ballast breakage in the form of corner abrasion and surface attrition is proposed. A series of DEM simulations are carried out on the large-scale cubical chamber to investigate the cyclic responses of ballast in the presence of fouling materials. The deformation and degradation of ballast during cyclic loading are investigated and reported from macroscopic and microscopic perspective. Detailed analysis on the micromechanical behaviour of fouled ballast in terms of inter-particle contact forces, inter-particle contacts within ballast depth and particle movement is provided, by which the impacts of fouling materials on ballast are explored and revealed.

Large-scale cubical box tests on fouled ballast
To investigate the effects of coal fines on the deformation and degradation of ballast, a series of large-scale cubical box tests was conducted by Indraratna et al. [14] using the Track Process Simulation Apparatus (TPSA). Figure 1a illustrates the schematic plane view and cross section of the TPSA, which can accommodate a ballast assembly 800mm long 9 600-mm wide 9 600-mm high. The wooden sleeper (200 mm (wide) 9 150 mm (high)) was supported by a 300-mm-thick layer of ballast assemblies, which was compacted to a unit weight of 15.3 kN/m 3 . Below the ballast layer was a 100-mm-thick layer of subballast and a 50-mm-thick layer of subgrade composing of coarse sand. Coal fines quantified by the Void Contamination Index (VCI) [50] were mixed with fresh ballast to represent a specific fouling level. The VCI is given by where e b and e f are the void ratio of ballast and fouling, respectively; G b and G f are the specific gravity of ballast and fouling, respectively; M b and M f are the dry mass of ballast and fouling, respectively. It is noted that coal fines used in the laboratory tests were chosen to be relatively dry owing to the prolonged exposure to dry weather in Australia [14]. Therefore, the cohesion and moisture content effects of fouling materials considered were not considered in their study. Cyclic loads were applied by a servo hydraulic actuator at a frequency of 15 Hz with a minimum pressure (q min ) of 45 kPa and a maximum pressure (q max ) of 420 kPa. The vertical settlement of the sleeper was recorded by an LVDT connected to the load actuator. The four vertical walls were allowed to move so that the predetermined lateral stresses could be maintained during the tests. Detailed information of the large-scale cubical box tests on fouled ballast can be found in Indraratna et al. [14].

Discrete element modelling of fouled ballast
In the present study, a DEM model was developed to simulate the cubical shear box test as reported by Indraratna et al. [14] using the Particle Flow Code 3D (PFC3D 5.0). A particle degradation subroutine was developed to capture surface abrasion and corner breakage of angular ballast in DEM. A series of cyclic loadings were applied to ballast assemblies to investigate the deformation and degradation behaviour of railway ballast at various fouling levels.

Model configuration
In the present study, a cubical box with dimensions of 600 mm (long) 9 300 mm (wide) 9 450 mm (high) was created as a representative section of ballast-sleeper as shown in Fig. 1b. In rail track, the dynamic loads induced by trains are transmitted through rails and fasteners to sleepers at a spacing of about 600 mm, which are embedded into ballast assemblies to uniformly distribute the concentrated loads from the steel rails to the ballast. It has been reported that the length of testing chamber (transverse) is expected to be greater than l/3 (l is sleeper length which is normally 2400 mm in standard gauge Australian tracks) considering the influence zone along the transverse direction (i.e. parallel to the sleeper) by exploiting double symmetry [41]. However, the current study focuses on the behaviour of ballast along the depth  and along the longitudinal direction (i.e. perpendicular to the sleeper). Therefore, only a section of 300 mm along the transverse direction was simulated in the current DEM model considering computational efficiency. The sleeper was modelled by a wall element (rigid element), and this approach has been used in previous studies [25,23]. A 300-mm-thick ballast layer was placed underneath the sleeper, having a particle size distribution that is identical to the one carried out in the cubical shear box test by Indraratna et al. [14]. In the DEM, the irregular shape of ballast has been found to be a critical factor influencing the mechanical performance and is directly associated with the inter-particle contacts and contact force chains developed during cyclic loading [33,53]. The approach to simulate angular-shaped ballast particles used in this study was adopted from Chen et al. [6], where the irregular particles were created using clumps logic with an image-aided DEM approach. It should be noted that the sphericity of ballast particles [54] in real track fields ranges from 0.65 to 0.95 [6]. In this study, the ballast particles are selected using sphericity as shape indicator to ensure that the particles represent the key shapes of ballast from each sieve sizes. Five ballast clumps as illustrated in Fig. 1c were generated inside the box at random particle sizes and orientations and then cycled to reach the initial porosity of 40%. The particle size distribution (PSD) of the ballast aggregate used in the current DEM model is shown in Fig. 2a [5,14,38,41]. For both clean and fouled assemblies, the total number of ballast aggregates was about 1235. Coal fines were simulated as spheres having sizes ranging from 0.1 mm to 5.0 mm (Fig. 2a). Four groups of fouled ballast assemblies were prepared at various fouling levels (i.e. VCI = 0%, 10%, 20% and 40%). The number of coal fines used in the simulations for VCI of 10%, 20% and 40% was 11,224, 22,550 and 45,147, respectively. A VCI of 40% indicates that 40% of the total ballast voids is occupied by the fouling materials as described by Tennakoon et al. [51]; it corresponds to a percentage fouling by mass of about 8% for coal-fouled ballast (the value is 2.4% and 4.5% by mass for VCI of 10% and 20%, respectively), whereby the granular mass is in a state of significantly impeded drainage requiring urgent maintenance. During simulation, the two vertical walls (wall 1 and wall 2) (Fig. 1b) were fixed to represent the plane-strain condition in the field (i.e. in longitudinal direction of track), while the walls parallel to sleepers (wall 3 and wall 4) were allowed to move so that a constant lateral confinement could be maintained. A lateral confinement, r y = 30 kPa, was selected to represent actual confining pressure provided by the shoulder ballast [24]. Cyclic loading (Fig. 2b) with a minimum value of q min = 45 kPa representing the in-situ stress state was applied onto the sleeper at a frequency of f = 15 Hz. The maximum applied load (q max ) varied from 230 to 420 kPa, which is typical stress levels that are induced by trains with axle loads of 25 t to 40 t in the field [49]. The limitations in the DEM (i.e. limited number of large asperities of ballast clumps and the overlooking of fatigue effect) have resulted in accelerating the simulation to stabilize quickly compared to the laboratory experiments [33]. Therefore, all the DEM simulations were carried out for 1000 cycles in the present study to capture the deformation and degradation behaviour of ballast during cyclic loading. The stress levels of lateral confinement and applied deviatoric load were achieved using a servo-controlled mechanism written in FISH language and developed by the authors specifically for the large-scale cubical box tests.

Ballast degradation model
In PFC, clumps are clusters of spherical particles that behave as an integrated granular assembly where their internal contacts (contact among particles within a clump) are being ignored during the numerical solution process, thereby the breakage of clumps is not allowed. To facilitate particle breakage, many alternative methods have been proposed, among which the most prevailing approaches are the Fragment Replacement Method (FRM) and the Bonded Particle Method (BPM). In the FRM, particle fragmentation is achieved by replacing the original particles with several self-similar sub-particles [20,31]. The main disadvantage of this method lies in the difficulties in selecting suitable particle breakage criteria and reasonable fragment replacement modes. In BPM, particle grains are treated as a group of sub-spheres that are connected by certain bonds where the failure of those bonds denoting particle breakage. This method has been successfully applied to the study of degradation behaviour of granular materials [37,42,52]. However, the highly irregular shapes of ballast particles require relatively small-sized spheres bonded together to represent the high angularity of particles, thus increasing the computational memory required for contact detection and solving force-displacement equations. These insufficiencies prevent the application of the BPM in studies where a large number of particles and long simulation terms are involved. In addition to the DEM-based particle breakage models, there are also some approaches attempting to combine the DEM with other numerical methods such as the FEM [39]. These coupled models (e.g. DEM-FEM [21,35], DEM-SBFEM [34] and DEM-XFEM [47]) have captured the particle breakage under different loading conditions. The basic theory of these combined methods is to discretize the polygonal particles (DEM) into a group of finite elements, with the stresses and the strains of these elements being calculated within the framework of the FEM to evaluate the breakage criteria. The polygonal particle assembly was then replaced with sub-polygons in the DEM to simulate the breakage of particle assembly. However, the breakage pattern of the simulated particle assembly is normally the bulk fracture, such as the rockfills in the Brazilian test. Ma et al. [36] explored the asperity yielding as well as the bulk fractures of angular rock grains under single grain crushing test by using combined DEM-FEM; however, unacceptably high computational memories are required for the simulations on ballast assemblies.
To capture the abrasion and breakage behaviour of ballast particles having irregular shapes, Liu et al. [30] proposed a two-dimension particle degradation model using PFC2D. In this study, the proposed model is further extended into three-dimension and implemented in PFC3D using FISH subroutines. The breakage criteria are based on the maximum contact force, and it is adopted to determine particle breakage in the existing DEM simulations [31]. The degradation (breakage) of a ballast particle can be categorized into body splitting and surface abrasion. Sun et al. [49] showed that the bulk fracture of ballast particles mostly occurs at extremely high loading frequency (f [ 30 Hz), while the predominant degradation mechanism of railway ballast is surface attrition and corner breakage. The developed subroutine used in the current simulation can be divided into three main modules as given in Fig. 3, including: (1) Potential abrasion pebble detection; (2) Contact force assessment; and (3) Abrasion pebble deletion and regeneration. The objective of the Module 1 is to detect all the potential abrasion pebbles which generally locate at the outer contours of ballast clumps. A radius threshold that controls the sizes of potential pebbles is adopted to avoid the abrasion of corners with large sizes and small curvatures. It is to be noted that the radius threshold ought to be selected with caution, as too low or too high a threshold value can lead to comprised or ill detection of corner pebbles and cause simulation inaccuracy. In the present study, the threshold is set to be 0.5 times the largest pebble radius in each clump by a series of preliminary trials. The size-dependent maximum contact force failure criteria are then incorporated into the Module 2, by which the degradation states of potential pebbles are determined. The allowable maximum contact force F c [27,30] for a potential pebble with diameter of d is specified as given by where the critical characteristic strength r c and the critical pebble diameter d c are used for strength normalization and particle size normalization, respectively; m is the Weibull modulus; a is the size scaling factor; P is the survival probability of the pebble and is randomly selected from a uniform distribution from 0 to 1. According to the analysis on particle shapes, the sizes of the potential pebbles located on the corners can basically represent the level of particle angularity. Therefore, the sizes of pebbles instead of the sizes of the irregularly shaped clumps are adopted in the breakage criteria. Pebbles whose maximum contact forces exceeded its F c are labelled as ready-to-break for the subsequent stage. The Module 3 is then triggered, where all the ready-to-break pebbles are automatically deleted from their parent clumps during the loading process, and newly created spheres are introduced to represent the particle fragments. Notably, the fragment spheres are in small sizes at their initial introduction to avoid the sudden change in contact forces with their neighbouring particles. The radii of newly created fragment spheres (R i ) depend on the non-overlapped volumes between the broken pebbles and their parent clumps (Vol i ), and the relationship is given by where fac s is the factor used for particle size scaling and was set as 0.2 in the current study by trial and error. The fragment spheres then linearly expand with a size expansion factor (fac e ) during the next 50 solving sequences to retain the mass conservation of whole system, as given by where N is the number of cycles for size expansion; R f is the ultimate radius of fragment spheres after expansion. By rewriting Eqs. (3)-(5), the fac e for each fragment sphere can be determined as given by It should be noted that the volume expansion might induce extra forces into the systems. Therefore, two methods were taken in the proposed degradation model to eliminate the induced forces: (1) the fragment spheres are created tangent to their parent clumps; this process can avoid the sudden change in contact forces the instant the creation; (2) the expansion cycles N need to be carefully selected. A too small value N accelerates the expansion process, but the created spheres will expand significantly during each expansion cycle and leads to extremely large contact forces with the neighbouring pebbles. However, a too large value of N could slow down the expansion process and the mass of the entire system could not be conserved during a short period. Detailed algorithm of the degradation subroutine can be found in Liu et al. [30].

Calibration of the proposed DEM model with direct shear test
To calibrate the DEM model, a series of large-scale direct shear tests simulations were carried out for ballast and coal fine mixtures. The dimensions of direct shear tests models were 300 mm (long) 9 300 mm (wide) 9 200 mm (high) similar to Indraratna et al. [15] (Fig. 4). Principally, parameters used to simulate ballast in DEM were selected by calibrating the DEM-based shear stress-strain responses with those measured from laboratory tests. These parameters were then varied interactively until the predicted shear stress-strain responses corroborate with those measured experimentally. In this study, a linear contact model [19] was adopted for contact resolution. The normal and shear stiffnesses for ballast (k nb ¼ k sb ¼ 1:7 Â 10 8 N/m) and coal fines (k nc ¼ k nc ¼ 2:8 Â 10 7 N/m) were used in the current DEM simulations [7]. It should be noted that the stiffnesses of particles from the same group (i.e. ballast or coal fines) were set the same irrespective of their sizes. The friction coefficients for ballast (l b ) and coal fines (l c ) and parameters for ballast degradation (a and r c ) were determined by comparing the simulation results with the shear stress-strain responses of fresh and 20%-VCI coalfouled ballasts in the laboratory tests [15]. Table 1 shows the detailed information of DEM simulations carried out for determining the key parameters (i.e. l b , l c , a and r c ). DEM-1 to DEM-8 were conducted for fresh ballast under a normal stress of 27 kPa, while DEM-9 to DEM-11 were for 20%-VCI assemblies. It can be seen from Fig. 4 that the selection of these micromechanical parameters has a significant influence on the macroscopic performance of ballast assemblies, with a larger shear stress of the assemblies being obtained at larger l b , a, r c and larger l c . It is seen that DEM-3 (l b ¼ 0:7, a ¼ 0:5 and r c ¼ 10 MPa) exhibits the best match with the laboratory experiments compared to other simulations on fresh ballast. Regarding the three simulations (DEM-9 to DEM-11) carried for 20%-VCI coal-fouled ballast at different values of l c , the stressstrain response obtained by DEM-10 (l c ¼ 0:3) was  Fig. 3 Diagrams of subroutine modulus for ballast abrasion slightly better despite there still exist some discrepancies as shown in Fig. 4d. It should be noted that parameter calibration is a cautious but cumbersome process, during which the mechanical responses on all fronts need to be comprehensively considered. Moreover, the effect of any individual micro-scale parameter on the mechanical response of the entire assembly is limited. Therefore, the selection of l c ¼ 0:3 is within the acceptable accuracy of the current DEM simulations. In addition, despite that the cubical box testing was carried out under cyclic loading conditions, particle displacements and velocities during simulations were small by comparing with the dimensions of the cubical shear box test. The complete loading process can be considered as quasi-static. In this regard, a local damping efficient, damp (a value of 0.7), was adopted for energy dissipation, and the introduction of this local damping can avoid the oscillation of small-sized particles so that the stability of whole system can be maintained.
The derived micromechanical parameters used in the current study are summarized in Table 2.

Cubical shear box simulation
Settlement of sleeper is an imperative factor affecting the stability of ballasted railway tracks. Figure 5a  times the total applied cycles for all the tested specimens, after which the settlements increase at an insignificant rate with number of cycles. The sleeper settlement for clean ballast when the aggregates were not allowed to break is also included in Fig. 5a. It is seen that the predicted settlement for unbreakable ballast assemblies shows lower value compared to the assemblies incorporating breakage. This observation indicates that the breakage of ballast exacerbates sleeper settlement, and neglecting the effect of particle breakage could lead to the underestimation of sleeper settlement in track. The ultimate sleeper settlement (S u ) for fresh and fouled ballast assemblies under different loading conditions is shown in Fig. 5b. In addition, results reported by previous studies for both fresh and coal-fouled ballast assemblies are also included in Fig. 5b to validate the current DEM model. It is shown that the predicted sleeper settlements are generally within the acceptable range compared to the reported laboratory test data. The S u is found to be linearly increasing with VCI. The increase in settlement with the increased fouling level has a detrimental effect on overall track stability [14]. In addition, the amplitude of loading (q max ) is found to have significant influence on the S u where the S u increases with the q max irrespective of the level of fouling. Subjected to cyclic loadings, the sharp corners of ballast break and the rough surfaces of particles abrade. Figure 6a presents the morphology evolutions of two typical simulated ballast particles captured at varied number of load  There are two predominant forms of particle breakage, i.e. surface attrition and corner breakage, have been captured in the present study. Notably, ballast particles become less angular with the increase in loading cycles (N) as shown in Fig. 6a. The decreasing angularity of the particle would decrease particle-to-particle interlock and therefore significantly reduces the bearing capacity of a ballast assembly. Figure 6b shows the predicted number of ballast fragments (B f ) with loading cycle (N) for assemblies at different fouling levels. It can be found that most particle breakage occurs within the 0.4 times the applied cycles, after which the assemblies experience slight increases in B f . This observation supports the fact that breakage of particles is a critical factor in increasing sleeper settlement when the assembly is subjected to repeated loadings. Figure 6c shows that the B f decreases with the increase in VCI as shown in Fig. 6c. The reduced ballast breakage for fouled assemblies is due to the 'coating effect' provided by the added fouling fines, which protect ballasts from abrasion during their rearrangement in response to cyclic loading [42]. In the 2D DEM simulations conducted by Indraratna et al. [14], the breakage of ballast particles was denoted by the number of broken bonds, which was similar to the number of fragments B f . It can be seen that the two indices show close relationships and can represent the ballast breakage indexes (BBI) from laboratory tests [14], which also validates the proposed DEM model in the present study.

Resilient modulus of ballast
During simulation, the applied cyclic stress (r cyc ) and the axial strain of sleeper (e a ) were monitored. The e a was determined by dividing the sleeper settlement by the depth of ballast layer, which was 300 mm in the present study. Figure 7a shows the predicted cyclic stress-axial strain responses of a fouled ballast assembly (VCI = 40%) recorded in the first five cycles for q max ¼ 420 kPa. It is   seen that the ballast assembly tends to respond elastically with increased loading cycles as reflected by the reduction in the size of the hysteresis loops, as also reported by past laboratory tests [42,3]. There exists irreversible axial strain during loading/unloading and the permanent axial strain decreases with the increase in load cycles, and the irrecoverable deformation is probably associated with particle rearrangement and breakage. The resilient modulus M R is an important parameter quantifying the deformation characteristics of ballast assembly. It is defined as the ratio of deviatoric stress (r cyc ) to the recoverable strain (e r ) during a load cycle. Figure 7b shows the 3D plot of the predicted M R and the fouling level (VCI) at different maximum cyclic stresses (q max ). As expected, M R increases slightly with the increase in q max for both fresh and fouled ballast assemblies. A significant drop in M R is observed when ballast is at higher fouling level, which implies the profound influence of the clogged fouling fines and the degradation of ballast particles on the resilient modulus of the granular assembly. Based on the simulated results, a relationship between M R , VCI and q max can be introduced as where a 1 , a 2 , b and c are dimensionless parameters depending on the level of fouling content and the applied cyclic loading conditions.

Degradation of ballast aggregates
In this study, the abrasion level is introduced to quantify the hierarchical process of ballast breakage during cyclic loadings. The abrasion level is defined as the number of times a particle experiences breakage during loading. A higher abrasion level denotes a more severe degradation of the particles that can lead to a large bulk fracture. Figure 8a shows a snapshot of abraded ballast grains for various abrasion level of a fresh assemblies for q max = 420 kPa. It is noted that particles that remain intact in the assembly are not shown in Fig. 8a. Most ballast particles experience less than three levels of abrasion, while a few particles have higher abrasion levels and are mainly distributed   Figure 8b and c presents the number of ballast particles at each abrasion level for fresh assemblies while at different loading conditions, and that for fresh and fouled assemblies for q max = 420 kPa. Statistically, most ballast particles undergo only one major level of abrasion, followed by relatively smaller amounts at Level 2 and 3. A limited volume of ballast experiences higher level of abrasion, especially for ballast at q max of 230 kPa. This is because, ballast particles gradually rearrange upon repeated loading and relative shearing between particle occurs causing abrasion (i.e. \ Level 3). After reaching a threshold stable state, any subsequent loading would cause less rearrangement and promote particle breakage in the form of either large corner abrasion or splitting across the body of particles. Moreover, the number of ballast aggregates at the initial abrasion level (i.e. \ Level 3) of fouled assemblies is reduced at low q max and high fouling level (VCI), which also verifies the 'coating effect' of fouling fines on ballast degradation. Figure 9 shows the images of broken fragments captured for fresh and fouled ballast assemblies at the end of loading for q max = 420 kPa. It is seen that the ballast fragments are mainly distributed within a trapezoidal region underneath the sleeper, albeit a small amount is observed in the area of crib ballast. This is expected as the ballast underneath sleeper majorly develop higher contact force than crib ballast during cyclic loading. It is noted that the identified trapezoidal region is consistent with the sleeper/ballast contact pressure zone marked during ballast vibration modelling by Ahlbeck [1]. The core degradation region containing more than 80% of the total broken fragments in assemblies is identified below the sleeper and is shown in Fig. 9. The pervasion angle (h) which is defined as the inclined angle of the boundary of core degradation region can be easily obtained by a subroutine developed in the FISH program. It is seen that the pervasion angle is significantly affected by the addition of fouling fines in assemblies. Fresh ballast shows the smallest pervasion angle (h ¼ 6 ), and the angle increases with an increasing magnitude of VCI. The increasingly dispersed ballast breakage underneath sleeper in fouled conditions not only increases the maintenance costs but also leads to the instability of track. The breakage of ballast along the depth underneath sleeper has also been investigated in the present study. Figure 10a shows the number of ballast fragments (B f ) at different depths for fresh ballast at various cycle ratio for q max = 420 kPa. With the increase of loading cycles ballasts along the depth experience an increased B f , while the most significant increase in B f is observed close to the sleeper owing to the high vulnerability to breakage of particles. There are some exceptions in B f at the late loading stage as shown in Fig. 10a, where the B f of upper aggregate slightly decreases and leads to the increase in B f for the deeper aggregate, and this observation could be attributed to fragment sedimentation during loading (i.e. broken particles moving downward). The displacements of ballast fragments at the end of loading are presented in Fig. 10b, where the fragments are represented as arrows showing their displacement directions. It is seen that most fragments move downward upon repeated loadings, which explains the reduction in B f for the particles close to the sleeper. In practical implications, the free drainage of ballast assembly would be significantly impeded by these sedimented fragments, and this observation also proves the detrimental effect of ballast degradation in track. Figure 11a shows the variation in B f with depth for fresh ballast when subjected to different q max . As expected, a greater B f is observed when the assembly is subjected to higher q max , and the B f decreases significantly with the increase in depth. A critical degradation depth (z c ) within which accumulates more than 80% of ballast fragments can be determined. For the tested fresh assemblies, the z c is about 0.6 times the total depth of ballast layer below sleeper, indicating this elevated vulnerability of particles within this depth. When it comes to the fouled assemblies, it is shown in Fig. 11b that the z c decreases to about 0.3 time the total depth for VCI of 40%. The decreased z c implies the concentration of particle breakage immediately below the sleeper in the presence of fouling, which further verifies the 'coating effect' of fouling materials on ballast degradation.

Microscopic analysis
To investigate insightfully the micromechanical behaviour of fresh and fouled ballast, a series of microscopic analysis are carried out to analyse the inter-particle contacts and the movement of ballast particles. Figure 12 shows the variations in the averaged inter-particle contact force (f ave ) and the averaged particle velocity (V ave ) of fresh ballasts with depth. Upon wheel loading, the stress imparted by sleeper/ ballast contact pressure attenuates with depth in the ballast layer, which has been verified by analytical models for track support systems [1,48] as well as by actual field measurements [13] as shown in Fig. 12a. The results of f ave for fresh assembly obtained by the present numerical models are also included in the figure for comparison. It is seen that both the f ave and the V ave gradually attenuate with the increase in depth. Above the critical degradation depth z c , i.e. 0.6 times the total depth of ballast layer, observes a strong degree of interlocking between particles. In terms of particle velocity, ballast within the 0.4 times the total depth exhibits relatively large velocities, indicating a considerably active destabilizing zone. Below this active zone, V ave of ballast gradually becomes smaller and almost remains stationary. The attenuated f ave and V ave within the ballast layer avoid the occurrence of large force concentration, while significantly reducing the possibility of relative sliding and rolling between ballast particles, thus leading to a reduction in particle breakage. Figure 13a presents the evolutions of the coordination number C n , which quantifies the average number of contacts per particle, of fresh and fouled ballast with loading for q max = 420 kPa. It should be mentioned that the interparticle contacts represent the ballast-to-ballast contacts solely in the present study. The C n for all the assemblies increases with N=N max and reaches its maximum value for N=N max = 0.4, thereafter the C n remains almost constant with N=N max . This has a direct relationship on the settlement behaviour of ballast. The sleeper settlement has found to increase slowly after N=N max ¼ 0:4 as shown earlier in Fig. 5. This observation implies that particle breakage becomes the dominant attribution to sleeper settlement after the initial cyclic densification of ballast assembly. The C n decreases with the increase in VCI. This is due to that fouling agents fill ballast voids and lead to the reduction in inter-ballast-ballast contact. The variation in the f ave with depth for fresh and fouled assemblies is shown in Fig. 13b. An increase in VCI decreases the f ave of assembly, which proves the aforementioned 'coating effect' of fouling fines, and this observation is consistent with the existing research [7,15]. Figure 14a presents the contact force distributions of fresh and fouled ballast at the end of loading for q max ¼ 420 kPa (N ¼ 1000 cycles). The contact forces are represented as lines connecting the centroids of ballast clumps with thickness being scaled by their magnitudes. The (a) (b) Fig. 12 a Variation in averaged inter-particle contact force with depths and b variation in averaged particle velocity with depths for fresh ballast Acta Geotechnica (2022) 17:3977-3993 3989 applied load is mainly supported by ballast underneath the sleeper with the inter-particle contact forces in this zone exhibiting relatively larger values compared to those located in the cribs. The significant contact forces underneath the sleeper are transmitted within a trapezoidal region. It is also noteworthy that the trapezoidal distributions of inter-particle contact force nearly coincide with the distributions of ballast fragments as shown in Fig. 9, which explains the emergence of core degradation region underneath sleeper. The anisotropy tensor [44] is adopted to describe the orientation of contact force, as given by where f k is the unit contact force vector, N c is the total number of contacts considered in the assembly, respectively. The anisotropy (a f ) and the principal orientation (h f ) of inter-particle contact forces can be determined considering the eigenvalues and the eigenvectors of the tensor [40], as given by where f 1 and f 3 are the largest and the smallest eigenvalues of f ij . h f is the elevation of the major principal force vector. Figure 14b shows the corresponding spherical histograms of contact forces of bottom ballast for fresh and fouled assemblies. It is seen that inter-particle contacts between bottom ballasts are nearly vertically distributed with the h f of about 90°. The anisotropy of contact force distribution (i.e. a f ) is 1.61 for fresh ballast, while it decreases with the increase in VCI. As the clogged fouling fines can provide sufficient capillary connections between ballasts, the imparted load forces are more uniformly distributed over the entire assemblies, therefore leading to the reduction in a f . This phenomenon also provides dependent evidence to the increased pervasion angle of core degradation region in fouled assemblies as shown in Fig. 9.

Limitations
Although the DEM model incorporating the particle breakage has captured the degradation and deformation responses of ballast assemblies during cyclic loading, the future studies should overcome the following limitations: 1. In this study, only five clumps of different shapes were used to generate ballast assemblies which may not completely represent the diverse angularities (shapes) of ballast in the track. Therefore, future simulations should focus on generating more angular ballast particles to investigate the mechanical behaviour of ballast assemblies during cyclic loading. 2. The current DEM model employs a linear contact model and considers the contact stiffness of ballast particles regardless of particle size. Hence, future DEM studies on assemblies of diverse particle sizes and angularities should be modelled using nonlinear Hertz model.

Conclusions
This study presents a 3D DEM model that is able to capture the degradation of irregularly shaped ballast particles in the form of surface abrasion and corner breakage. The developed particle degradation model adopted the maximum contact force failure criteria considering the Weibull distribution effect and the size hardening law on ballast strength. DEM simulations were carried out for a cubical (a) (b) Fig. 13 Coordination numbers and the averaged inter-particle contact force f ave of fresh and fouled assemblies with loading at q max of 420 kPa shear chamber incorporating a representative track section by implementing a particle degradation subroutine developed in FISH language in PFC3D. The parameter calibration was conducted by corroborating the predicted shear stress-strain responses of ballast under direct shear conditions with the measured laboratory data. Cyclic loadings with different amplitudes were applied onto fresh and fouled ballast assemblies at a frequency of 15 Hz, and the deformation and degradation of ballast for various levels of fouling were captured. The main conclusions can be summarized as follows: 1. Significant ballast breakage (i.e. the number of frag-mentsB f ) accompanied with a rapid increase in sleeper settlement (S) is observed at low number of cycles; thereafter, the ballast exhibits a plastic shake-down phenomenon with a stably increasing sleeper settlement. At q max ¼ 420 kPa, fresh ballast experiences an ultimate settlement of the sleeper (S u ) and a B f of 19.4 mm and 200, respectively; the S u increases to 25.3 mm and the B f decreases to 103 as VCI increases to 40%. The S u and the B f of fresh and fouled ballast increase with the increase in the applied loading amplitudes (q max ). 2. The resilient modulus (M R ) is directly affected by particle breakage and fouling fines. At q max ¼ 420 kPa, the M R is 151 MPa for fresh ballast, and it decreases to 85 MPa when VCI increases to 40%, and the values decreases to 132 MPa and 75 Mpa, respectively, with q max decreasing to 230 kPa.
3. The study proved that much of the ballast experiences less than three encounters of abrasion in response to their initial rearrangement, after which little breakage (abrasion level [ 3) is induced by the subsequent loadings. Due to the 'coating effect' provided by the fouling coal fines, the coordination number (C n ) and the averaged inter-particle contact force (f ave ) of the assembly reduce from 5.8 to 5.0, and 8.3 to 5.2 kN (measured underneath the sleeper at q max ¼ 420 kPa) when VCI increases from 0 to 40%, and this protects excessive breakage of ballast. 4. More than 80% of ballast fragments are mainly distributed within a trapezoidal region underneath sleeper. As VCI increases from 0 to 40%, the region becomes increasingly dispersed in assemblies with the pervasion angle (h) increasing from 6°to 22°due to the more uniformly distributed inter-particle contact forces (the anisotropy coefficient a f decreases from 1.61 to 1.45). 5. The breakage of ballast mitigates with the increase in depth below sleeper due to the attenuated inter-particle contact force and particle moving velocity. Beyond a depth of 180 mm, the damage to fresh ballast is insignificant, and this critical depth decreases to around 90 mm for ballast at VCI of 40%. Upon repeated loading, ballast aggregates at all depths experience an increased breakage accompanied by sediments from broken fragments.  Fig. 14 a Distribution of inter-particle contacts between ballast in fresh and fouled assemblies at q max of 420 kPa and b spherical histograms of inter-particle contact distributions for bottom ballast Acta Geotechnica (2022) 17:3977-3993 3991 Funding Open Access funding enabled and organized by CAUL and its Member Institutions.

Data availability
The data used to support the findings of this study are available from the corresponding author upon request.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.