1 Introduction

It is commonly accepted that, compared with railway ballasted tracks, well-constructed slab-tracks generally result in significantly lower maintenance needs for restoration of track geometry due to settlements. However, for settlements over a certain threshold, possibilities for routine un-expensive maintenance are much more limited. For granular embankments, a certain level of uncertainty on the development of settlements due to traffic loading exists within the embankment fill itself. This is due to the complexity of granular matter, mainly stemming from its discrete nature. This discrete nature becomes especially relevant when particle sizes become meaningful in relation to the considered structure. This is the case with rockfill embankments, i.e. embankments constructed with blasted and/or crushed rock of relatively large particle size. In this regard, the so-called discrete element method (DEM) allows to explicitly consider the discrete nature of the material. The effect of processes at particle level on the macroscopic response can thus be investigated and insight into the principles of material behaviour provided.

There is a lack of studies regarding full-scale measurements of settlements for rockfill railway embankments. Empirical studies involving railway embankments are generally performed on reduced scale models, e.g. [1, 2]. Full-scale studies on high rockfill embankments are mainly for dams, e.g. [3, 4]. Although these infrastructures present fundamental differences with railway embankments, dam height was long identified as one of the most influencing factors defining settlements [5]. Nevertheless, apart from this lack of specific full-scale measurements, empirical studies alone fail to provide a deeper insight into the fundamentals of granular matter due to limitations in conducting comprehensive measurements at particle level. DEM grants extensive access to information at particle and contact level in an unprecedented way, allowing to study and explicitly and quantitatively characterize granular assemblies in a way that cannot be matched by physical tests, as for example in connection with soil fabric and force chains [6,7,8,9].

Several modelling examples of railway ballast using DEM can be found. There are studies showing the importance of particle shape in obtaining a more realistic response, like in [10, 11]. In particular, these studies showed the need for clumps of spheres as opposed to individual spheres in order to capture the effect of interlocking on rearrangement. Studies also show that simple clumps can result in realistic results [11,12,13] presenting a very computationally efficient alternative to complex irregular shapes trying to reproduce more realistic shapes for ballast particles. The importance of particle degradation when modelling ballast under triaxial loading has also been investigated [11, 12, 14]. There are even examples where the ballast layer has been modelled showing the significance of particle degradation and fouling on settlements [15, 16]. It should be observed that all the above numerical studies refer to the ballast layer, typically 30–40 cm thick measured from the underside of the sleeper. For the case of embankments, the fill between the track structure and the foundation is usually measured in meters, with values up to 10 m and above not being uncommon.

Despite all the above advances in numerical modelling, to the best of the authors’ knowledge, no attempt to model a full railway embankment under traffic loading, that besides includes particle degradation, has been undertaken to date. This is mainly due to the great computational time demands associated with modelling large collections of particles or large numbers of load cycles using DEM. This becomes aggravated by including degradation, especially if using bonded spheres as in some of the above studies. Moreover, bonds are governed by numerous parameters lacking a direct connection with physically-measurable material properties. Due to the mentioned lack of full-scale studies on rockfill embankments, these parameters need to be entirely determined by calibration under conditions different to the ones intended for, e.g. triaxial conditions. This renders the model rather vulnerable outside those conditions, i.e. for the actual modelling of an embankment.

Recently, a DEM material model intended for modelling large constructions has been presented [17]. Particles consisted of simple breakable clumps of spheres including asperities. Both corner breakage and particle splitting were implemented and differentiated depending on local contact geometry. Moreover, these processes were controlled by fewer parameters and with a clearer physical meaning than in the case of bonded spheres. All in all, it resulted in more computationally efficient and robust models ideal for their application for high rockfill embankments. This article presents the implementation of the breakable clumps in the actual modelling of rockfill railway embankments. Figure 1 shows a schematic representation of a possible embankment cross-section for slab-track configuration. Only the coarsest fraction of the rockfill has been targeted to reduce computational time. In the model, embankments are generated by successive dumping and compaction of layers of clumps on top of each other mimicking the construction of real embankments. After generation, permanent loading representing the weight of the infrastructure above the unbound layer (see Fig. 1) is applied followed by activation of particle degradation and cyclic loading representing traffic on top of the permanent load. The model is used to study the influence of embankment height and particle breakage on the mechanical response to cyclic loading, especially on settlements accumulation and resilient response. The main unique contribution of this article is the use of DEM for modelling high railway rockfill embankments under traffic loading, furthermore including degradation, something not attempted to date.

Fig. 1
figure 1

Schematic cross-section of modelled embankment for slab-track configuration. Dimensions in mm

Full size image

2 Material model

The main features of the breakable clumps presented in de Frias Lopez et al. [17] are briefly summarized here. Particles consisted of simple breakable clumps of spheres formed by a tetrahedral body with four asperities (see Fig. 2). Breakage was implemented by releasing selected spheres from the clump. Spheres representing asperities were released to simulate the breakage of angular projections, i.e. corner breakage; all spheres pertaining to a clump were released to represent particle splitting. In particular, the latter resulted in two distinctive fragment sizes (see Table 1), a type of splitting observed in real ballast particles [16, 18]. Breakage processes were controlled by the values of the acting normal contact force \({f}_{\mathrm{n}}\) and the particle loading configuration, characterized by the clump coordination number \({C}_{\mathrm{n}}\), i.e. number of contacts with neighbouring particles. For every calculation step, \({f}_{\mathrm{n}}\) values were compared with the value of the maximal recorded normal contact forces \({f}_{\mathrm{n},\mathrm{max}}\) (generally set to zero on the onset of loading) for each of the contacting entities. If \({f}_{\mathrm{n}}>{f}_{\mathrm{n},\mathrm{max}}\), the degradation algorithm activated distinguishing between the cases presented in Table 1, where \({f}_{\mathrm{n},\mathrm{u}}\) and \({C}_{\mathrm{n},\mathrm{min}}\) are the parameters defining the triggering of breakage. The former parameter represents the value of loads liable to produce fracture, whereas the latter represents a minimum value of \({C}_{\mathrm{n}}\) above which clumps are no longer liable to experience splitting. In de Frias Lopez et al. [17], \({f}_{\mathrm{n},\mathrm{u}}\) values were assigned following a log-normal distribution representative of the experimental probabilistic distributions for indentation testing of granite reported in [19]. The value of \({C}_{\mathrm{n},\mathrm{min}}\) was set to 3 representing a configuration equivalent to an indirect tensile strength tests as in [16, 18]. Monotonic and cyclic triaxial tests were performed under a range of stress levels. A comparable degree of resemblance to empirical results as previous state of the art modelling efforts with bonded asperities was observed.

Fig. 2
figure 2

Regular tetrahedral clump with four asperities

Full size image
Table 1 Summary of degradation algorithm
Full size table

It should be noted that the material model presented in [17] is in line with recent approaches advocating for simpler clumps where simplifications are compensated, to some extent, by contact laws where damage processes are included [13]. Furthermore, the model was validated using different tests types and multiple measurements, i.e. monotonic and cyclic triaxial testing using results for axial and volumetric deformations as well as resilient properties. All in all, it resulted in an efficient and robust material model liable to be implemented in the present modelling of rockfill embankments.

3 DEM modelling

Modelling has been undertaken using the PFC3D 4.00 DEM software environment [20]. Particles were represented by simple tetrahedral clumps with four asperities (Fig. 2) following [17]. Half-track embankment fill models were generated for different embankment heights and cyclically loaded to simulate the effect of traffic. Results at both macroscopic and particle level were obtained and analysed.

3.1 Embankment generation

When modelling three-dimensional ballasted tracks using DEM, it is common practice to only consider half-tack [15, 21]. Regarding the length of the model in the track direction, there are examples of one-sleeper [15], three-sleeper [21] and even of up to six-sleeper [22] long models. In the present study, a four-sleeper half-embankment model was chosen, equivalent to 2.6 m along the track (650 mm sleeper spacing assumed). Cross-section dimensions were based on Fig. 1 for target heights \(H\) ranging from 2 to 10 m. The term “target” was used as, due to the implemented generation procedure (see below), the exact value of the dimensions for the compacted embankment could not be predicted beforehand.

Particle size and other micromechanical properties were generally chosen following [17] where parameter values had mostly been defined based on previous numerical and experimental studies [11, 12, 15, 23, 24]. All vertical boundary symmetry planes were modelled by planes with same stiffness and friction as the initial average particle values. This represented the effect of the surrounding material on the modelled portion of the embankment (see Sect. 3.2). Same applied to the bottom wall representing the foundation. The main input material properties are summarized in Table 2. A local damping coefficient of 0.7 was applied throughout the simulations following [17] in order to absorb vibration energy and reach mechanical equilibrium within a reasonable number of computational cycles (see Sect. 3.2).

Table 2 Main input micromechanical parameters for embankment generation
Full size table

Embankments were generated by successive dumping and compaction of layers of unbound material on top of each other until obtaining a compacted height of at least the chosen target height. First, a rectangular parallelepiped container of six walls was generated. These walls represented the overall boundaries of the embankment, including symmetry planes and the unbound layer/subgrade interface. An example of overall material vessel is shown in Fig. 3. An auxiliary trapezoidal container for generation and dumping of the first layer was generated within the overall container. The dumping container made use of the front and back boundary walls together with the half-track symmetry plane, so only three walls were created: a top and bottom horizontal walls together with an outer lateral wall with same slope as the intended embankment lateral slope, i.e. 1V:1.5H (see Fig. 4). The thickness of this container \({t}_{\mathrm{cont}}^{\mathrm{p}}\) was set as the approximately expected thickness of the layer after compaction. A value of 0.5 m, which can be representative of rockfill embankments, was chosen. For the first layer, the bottom width of the dumping container was set as the embankment targeted bottom width for the chosen \(H\) following Fig. 1 and correction for the layer width based on the embankment slope and the compacted height of the so far generated embankment was made for successive layers. The dumping container was first generated with its bottom wall at a distance of \({h}_{\mathrm{cont}}^{\mathrm{p}}\) from the bottom of the overall container. A value of 0.1 m was chosen. Same distance \({h}_{\mathrm{cont}}^{\mathrm{p}}\) was kept from the top of the last compacted layer to the bottom of the dumping container for successive layers. The dumping container was filled with a collection of randomly placed spheres according to a random-number generator seed with a porosity \({n}_{\mathrm{cont}}^{\mathrm{p}}\) following the particle size distribution defined in Table 2. The porosity was set to 0.5 to avoid high initial contact forces due to overlapping. Ball and wall properties were assigned with balls kept as frictionless as the mixture is relaxed. Balls were then substituted by frictionless clumps and the mixture was relaxed to equilibrium. Particle friction was applied and the three additional walls defining the dumping container were deleted to allow for free fall under a gravity acceleration of 9.81 \(\mathrm{m}/{\mathrm{s}}^{2}\) for a total of 1.0 s, allowing all material to settle under gravity on top of the so far compacted embankment.

Fig. 3
figure 3

3D view of final compacted embankment for \(H=3\) m about to be loaded by a loading wall (see Sect. 3.2). Different colours indicate different clump layers

Full size image
Fig. 4
figure 4

Embankment generation of layer 4 for \(H=2\) m. From above to below: compaction of layer 3 by a compacting wall; layer 4 generated in a container above layer 3 and about to be dumped; layer 4 dumped and settled over layer 3 and about to be compacted. Different colours indicate different clump layers

Full size image

Once the dumped layer had settled and furthermore had been relaxed, a compacting wall with same friction and stiffness as the bottom wall was created. The compacting wall was created right on top of the highest clump and progressively accelerated downwards to a maximum compacting velocity \({v}_{\mathrm{max}}^{\mathrm{c}}\) of 0.01 m/s, a rate of loading low enough to ensure substantial rearrangement. Once a maximum normal stress \({\sigma }_{\mathrm{n}}^{\mathrm{c}}\) of 100 kPa was achieved, the wall was halted to a stop. This value can be considered as a representative vertical stress level near the layer surface in the proximity of a roller compactor [25, 26]. Unloading was implemented by reversing the direction of the compacting wall in a similar fashion as loading until the measured reaction force on the compacting wall was null. This loading–unloading process was repeated for a total of five compaction cycles \({N}_{\mathrm{cl}}^{\mathrm{c}}\). On completion of compaction the wall was halted and inactivated by setting its stiffness to zero. The purpose of the implemented compaction procedure was to result in a similar level of compaction and internal stresses as a compacted granular layer rather than to fully reproduce the mechanics of compaction of real granular layers. Table 3 summarizes the main input parameters defining the dumping and compaction processes whereas Fig. 3 shows an example of compacted embankment for \(H=3\) m. Figure 4 schematically shows the intermediate stages for layer generation by dumping and compaction.

Table 3 Parameters for material dumping and layer compaction during embankment generation
Full size table

3.2 Embankment loading

After generation, the embankment was loaded with a permanent load equivalent to the self-weight of the infrastructure above the granular layer, estimated as roughly 25 kPa based on the geometry shown in Fig. 1, together with unit weight values of 25 kN/m3 for bound layers and 23 kN/m3 for unbound ones. Permanent loading started by creating a horizontal loading wall on top of the compacted embankment (see Fig. 3) with a lateral width of 4.2 m, equal to the rockfill model top width (see Fig. 1), and same stiffness and friction as the average particles values. The wall was progressively accelerated downwards to a maximum velocity \({v}_{\mathrm{max}}^{\mathrm{PL}}\) of 0.01 m/s, a rate of loading equivalent to \({\dot{\varepsilon }}_{\mathrm{a}}=0.005\) s−1 for \(H=2\) m as in [17]. This ensures that results for heights \(H\ge 2\) m are load rate independent (more on this below). Once the wall normal stress reached 25 kPa, the wall was halted to a stop for 5000 calculation steps to ensure equilibrium. Table 4 summarizes the main input parameters defining the permanent loading procedure.

Table 4 Main permanent and traffic loading parameters for embankments
Full size table

Traffic load was based on a single axle load of 300 kN for heavier than normal rail traffic at low to moderate speeds, i.e. negligible dynamic effects at rockfill level. This load was distributed over two to three sleepers and further spread down to the interface between the frost protection layer and the rockfill level (Fig. 1), resulting in an equivalent surface similar to the loading wall surface with and acting pressure of 13.5 kPa. Consequently, the implemented axle load should not be regarded as uniformly distributed over the track length but rather over the loading wall length. Considering that the main focus of this article lays on the settlements due to traffic loading, friction is added to the vertical boundaries to account for the effect of the material surrounding the studied portion of the track and even include a certain shear reversal effect, known to affect settlements [27]. Nevertheless, how to handle boundaries at macroscopic level is not always clear beforehand when using DEM, as macroscopic symmetry does not necessarily translate into symmetry at the particle level in a discrete system. In particular, a flat boundary (vertical wall) differs from reality in this case, where interlocking with other particles at the assumed boundaries would take place. Furthermore, adding friction at the boundaries is observed as common practice in several studies on ballasted tracks or box tests [10, 22, 28, 29].

Traffic loading was applied on top of the permanent loading for a total of 500 load cycles (10 conditioning cycles followed by 490 full load cycles). The loading started with the loading wall being progressively accelerated downwards to a maximum velocity \({v}_{\mathrm{max}}^{\mathrm{TL}}\) of 0.01 m/s. Once a maximum total normal stress of 38.5 kPa had been achieved the wall was halted to a stop. Unloading was implemented by reversing the direction of the loading wall in similar fashion as loading until the measured reaction force on the compacting wall was 25 kPa. Table 4 summarizes the main input parameters defining the traffic loading.

Regarding the selected loading rate and damping coefficient, it has been empirically shown that, for conventional speeds, load frequency does not significantly affect the accumulation of permanent strains in ballast [30, 31]. If this is true for the ballast layer, i.e. granular layer directly under rolling traffic, it can be considered as even more relevant for the rockfill, i.e. at a certain depth from the track. Dynamic effects are expected to be further attenuated by intermediate layers, like frost protection layers, and hence of a lesser magnitude. This, together with the extensive computational time demands when modelling true dynamic processes using DEM, leads to perform simulations as load rate independent, i.e. slow low rate. Under such conditions, the real dissipative mechanisms need not be replicated in DEM, therefore the model can be damped using the default local damping coefficient of 0.7 [20]. This is deemed as a reasonable simplification that significantly reduces computational time, making feasible to work with large numbers of particles under large numbers of loading cycles.

A series of test not including degradation, i.e. “unbreakable”, and another series including degradation, i.e. “breakable”, were performed, where the degradation mechanism was activated on the onset of traffic loading. Following [17], \({f}_{\mathrm{n},\mathrm{u}}\) values for the breakable test series were assigned following the log-normal distribution in Fig. 5, representative of the reported results for indentation testing of granite in [19]. The value of \({C}_{\mathrm{n},\mathrm{min}}\) was set to 3 following previous studies [16,17,18].

Fig. 5
figure 5

Distribution of fracture loads \({f}_{\mathrm{n},\mathrm{u}}\) for breakable tests series

Full size image

4 Results and discussion

4.1 Compaction

Embankment generation for target heights \(H\) ranging from 2 to 10 m in steps of 1 m was carried out. The basic volumetric properties after generation are summarized in Table 5. Final compacted embankment heights \({H}_{\mathrm{emb}}\) were quite close to the target ones, as in general all other dimensions were. The number of clumps varied from less than 5000 to clumps in excess of 50,000, highlighting the computational time challenges associated with higher embankments. Figure 6 shows the particle size distribution for selected embankments together with the one for the numerical triaxial tests in [17]. All embankments follow almost identical nearly-uniform distributions, being the gradation curve shapes very similar to the ones for the numerical triaxial tests. Embankment porosities, computed as the average value over the entire volume of the embankment, are presented in Fig. 7. An average decrease of porosity with embankment height can be observed together with a slight level of oscillation. Lower porosities for higher embankments are indicative of a slightly more compacted state (see below). Average embankment porosities are not far from the 0.51 value after triaxial specimen generation in [17]. This is expected to be the case especially for the areas closer to the half-track symmetry plane and away from the free slope, where the porosity is expected to be lower than the average values in Fig. 7. Considering this and the similarities in particle size distributions, a reasonable level of similarity between the triaxial specimen in [17] and the embankment conditions far from the free slope can be established.

Table 5 Embankment volumetrics on generation completion
Full size table
Fig. 6
figure 6

Particle size distribution for compacted embankments. Triaxial results after numerical tests in [17]

Full size image
Fig. 7
figure 7

Porosities of compacted embankments

Full size image

Figure 8 shows the average horizontal normal stresses after embankment generation for all tests together with the theoretical at-rest values (see below). Numerical stresses are obtained as average values for the normal forces acting on the different vertical confining walls. The average for the front and back walls corresponds to the longitudinal confinement and the half-track symmetry wall to the lateral. Both lateral and longitudinal components show similar values with nonetheless slightly higher lateral stresses. Although regarding the longitudinal direction, clumps are confined between two walls, i.e. front and back walls, a slightly lower value of the resulting average confinement was not fully unexpected. Averaged longitudinal stresses partly account for the material close to the free slope region, where horizontal components caused by overburden or otherwise locked stresses resulting from compaction are expected to be lower due to free slope boundary and lateral spreading. For the lateral confinement, the free slope can be expected to have a lesser influence on horizontal pressures due to overburden at the half-track symmetry plane, as this is the plane furthest away possible from the free slope.

Fig. 8
figure 8

Embankment horizontal stress on completion of compaction. Numerical results for longitudinal and lateral components; theoretical at-rest horizontal earth-pressure values after Eq. (3)

Full size image

Confinement values seem to grow in a nearly-linear fashion with embankment height, which was anticipated as result of an expected linear growth of the overburden pressure and consequent growth of the horizontal component at rest as follows. It is common practice in soil mechanics to determine the ratio between the horizontal and vertical stresses at rest, i.e. the at-rest earth pressure coefficient \({K}_{\mathrm{o}}\), for granular soils by using the empirical expression [32]:

$${K}_{\mathrm{o}}=1-\mathrm{sin}(\phi )$$
(1)

where \(\phi \) refers to the drained friction angle. This angle can be obtained by performing monotonic triaxial tests to failure as in [17]. The results of the monotonic tests at failure for the breakable assemblies in [17] are summarized in Fig. 9 in the p-q space, i.e. mean normal stress p and deviatoric stress q plane. By means of linear regression of the form \(q=M\cdot p\), where M represents the slope of the failure line in the p-q space, a value of 1.74 is obtained for M (\({R}^{2}=0.9986\)). The relation between \(M\) and \(\phi \) is commonly expressed as:

$$M=\frac{6 \mathrm{sin}(\phi )}{3- \mathrm{sin}(\phi )}$$
(2)

where \(\phi \) can be backcalculated to an angle of approximately 43°. The average horizontal pressure can then be determined as:

$${\sigma }_{\mathrm{h}}={K}_{\mathrm{o}}\cdot {\sigma }_{\mathrm{v}}+{\Delta \sigma }_{\mathrm{h}}={K}_{\mathrm{o}}\cdot \frac{{H}_{\mathrm{emb}}}{2}\cdot {\gamma }_{\mathrm{s}}\cdot (1-{n}_{\mathrm{emb}})+{\Delta \sigma }_{\mathrm{h}}$$
(3)

where \({\sigma }_{\mathrm{v}}\) represents the average vertical component of the stress caused by the overburden due to self-weight of the material, which is obtained based on the particle unit weight \({\gamma }_{\mathrm{s}}=25\) kN/m3 (see Table 2) and the embankment porosities \({n}_{\mathrm{emb}}\) in Fig. 7. \({\Delta \sigma }_{\mathrm{h}}\) represents the overall additional effect of the locked-in stresses due to compaction as otherwise Eq. (1) would underestimate at-rest values for compacted granular materials, e.g. [33]. By roughly setting \({\Delta \sigma }_{\mathrm{h}}\) to 15 kN/m2 (this value is justified later in Sect. 4.2), a good agreement with the confinement results in Fig. 8 is obtained. Nevertheless, independently of the assumed value of \({\Delta \sigma }_{\mathrm{h}}\), a friction angle of 43° based on the triaxial tests is able to capture the overall trend in increasing confinement with embankment height.

Fig. 9
figure 9

Shear strength in the p-q space after monotonic triaxial tests to failure in [17]. Linear regression line \(q=M\cdot p\) included as dashed line

Full size image

4.2 Traffic loading

Settlements due to traffic load for both the breakable and unbreakable series are summarized in Fig. 10, where embankment height refers to total height after application of the permanent load \({H}_{\mathrm{emb}}^{\mathrm{PL}}\). Given that the embankment generation was based on dumping containers filled by randomly placed spheres according to a random-number generator seed (see Sect. 3.1), each embankment can be regarded as a stochastic realization for the considered target height \(H\) and controlled by the selected seed. As an example of this stochastic nature, the whole generation and loading procedures were repeated for selected embankments by changing the generator seed. These can be seen as new realizations of the considered target heights and are referred to as Seed 2 in Fig. 10. The newly generated embankments result in different settlements to the original embankments, displaying the aforementioned stochastic nature.

Fig. 10
figure 10

Settlement of the embankment due to traffic load for breakable (Brk.) and unbreakable (Unbrk.) tests

Full size image

Settlements in Fig. 10 are characterized by a noticeable degree of dispersion and seem to be unaffected by embankment height. A substantial effect of breakage is nevertheless generally observed, where results do not seem to follow any discernible pattern. However, for the majority of cases, settlements for breakable assemblies are higher than for unbreakable and in occasions even doubling the latter. This effect of breakage is especially noticeable for intermediate and high embankments, where dispersion is higher, i.e. increased uncertainty in developed settlements.

Given the stochastic nature of the settlements for any given embankment height, these can be studied by means of statistical test. In particular, a matched paired-sample t-test can be conducted to determine if breakage is statistically significant. Samples are paired by embankment height and generation seed, i.e. by generated embankment, resulting in 12 pairs of samples. The considered null hypothesis is that breakage has no effect on settlements. Results are summarized in Table 6. The null hypothesis is therefore rejected as a statistically significant effect of breakage is observed.

Table 6 Summary of results of paired-sample t-test to study the effect of breakage in settlements
Full size table

The evolution of settlements after conditioning with the number of loading cycles for selected embankments of breakable and unbreakable assemblies is shown in Fig. 11. Embankments do not reach shakedown, i.e. a decreasing incremental plastic strain per load cycle until eventually approaching zero, which would have resulted into settlements eventually reaching a flat response. This type of response deviates from the one observed in ballast materials subjected to low deviatoric to confinement stress ratios under triaxial conditions [24], later reproduced numerically in [12, 17], where shakedown was generally achieved. This indicates that the fundamentally different boundary conditions between the triaxial configuration and the free slope boundary in embankments play a significant role on macroscopic response under cyclic loading. For breakable assemblies, some embankments seem to follow a logarithmic growth, i.e. target heights \(H\) of 4 and 6 m, as long had been observed for granular materials [34, 35]. The rest seem to follow either exponential growth, also common in granular materials [36] and most clearly observed here for \(H=10\) m, or something in between exponential growth and two differentiated phases of logarithmic growth marked by an relative acceleration in settlements after approximately 150 load cycles. No clear influence of embankment height in dictating the type of behaviour is observed. For unbreakable tests, settlements seem to mostly follow exponential trends. Again, no clear influence of embankment height is observed. It can also be observed that initial values, i.e. settlements after completing the conditioning phase, are rather similar for the unbreakable tests and trends in their evolution are not so dramatically different between embankments. For breakable tests, a greater dispersion in initial values and trends is observed, in agreement with the already observed increased uncertainty in final settlements attributed to breakage in Fig. 10. It is remarkable that although breakage generally results in higher settlements (Fig. 10) it can also, in occasions, result in a more stable logarithmic behaviour with reduced settlements. Regarding the number of modelled loading cycles, de Frias Lopez et al. [17] suggested that rearrangement and degradation develop considerably faster in numerical specimens and hence a reduced number of load cycles could represent a larger one in reality (a ratio of approximately 1 to 100 was observed in their triaxial tests). This was suspected to be caused by the chosen loading rate (see Sect. 3.2) maximizing rearrangement together with the limited number of large asperities and ignoring the effects of fatigue and slow crack growth as earlier suggested by Lu and McDowell [12].

Fig. 11
figure 11

Evolution of embankment settlements after conditioning under traffic load for a breakable and b unbreakable tests

Full size image

The resilient modulus can be estimated based on the stress–strain response under cyclic loading. Moduli are estimated as the secant modulus during unloading, as shown in Fig. 12. The resilient response on completion of traffic loading, i.e. based on cycle 500, for all embankments is summarized in Fig. 13. A linear stiffening effect of embankment height can be observed for both the breakable and unbreakable assemblies. This is in agreement with results for triaxial tests in the sense of linear growth of the modulus with confinement [24]. Furthermore, it can be observed that degradation has a minor effect on stiffness, as was the case in [17]. It was speculated in the latter that, due to the irreversible nature of breakage, its effect is more profound in the permanent component of the strains than in the resilient, as seems to be the case in this study. It can also be observed that breakable assemblies are generally slightly stiffer than unbreakable ones. In fact, low levels of degradation could be expected to have a slightly compacting effect by facilitating rearrangement into more compact configurations. This could also explain why, for very few cases, settlements for breakable embankments are smaller than for unbreakable ones (Fig. 10). In such cases, it can be inferred that the beneficial stabilizing effect of breakage overcomes the additional settlements it produces. For the triaxial tests in de Frias Lopez et al. [17], breakage was considerably more extensive (see below) and therefore not necessarily resulting in more compacted specimens. In fact, a profound effect of breakage on the bearing skeleton of the specimens, characterized by the coordination number, was observed.

Fig. 12
figure 12

Stress–strain response due to traffic load at the loading wall after loading cycle 100 for breakable embankment \(H=5\) m, breakable assembly. Resilient modulus can be estimated as the secant modulus during unloading, i.e. the tangent for the dashed line

Full size image
Fig. 13
figure 13

Embankment resilient modulus after traffic loading for breakable and unbreakable tests

Full size image

Regarding results at particle level, Fig. 14a shows the final accumulated degradation for the breakable tests. Corner breakage and particle splitting are presented as number of events relative to the initial number of clumps. For all embankments, corner breakage is dominant compared to particle splitting, as was the case for the triaxial tests in [17]. Degradation values for embankments are clearly below the ones for the triaxial tests, where corner breakage in excess of 15–20% was generally observed. This can be explained by the considerably lower stress levels for embankments. No clear correlation between the evolution patterns of corner breakage in Fig. 14b and settlement development in Fig. 11 is found, contrary to what was the case for the triaxial tests in [12, 17]. This could be in principle attributed to relatively low levels of degradation for embankments. Therefore, the way degradation dictates patterns in the macroscopic response may not be as distinct as for the triaxial tests. In fact, degradation levels are so low that its impact on the gradation was nearly imperceptible. The only significant effect was a new additional finer fraction consisting of spheres resulting from corner breakage and particle splitting, thought this fraction represented less than 0,05% of the total volume for all cases. Finally, an average oscillatory increase in degradation with embankment height is generally observed, thought the highest values of degradation are for the lowest embankment. This could in principle be explained by boundary effects and singularities on the contact forces network that are overrepresented due to the lower number of clumps in contact with the walls. Higher than anticipated contact forces can develop for a number of clumps, where these become locked in an unfavourable position. This explanation is in agreement with previous studies showing the significance of individual particles trapped in singular configurations on the force transmission for contact networks of limited magnitude [9, 37].

Fig. 14
figure 14

Embankment degradation under traffic loading for breakable tests: a effect of height on final values of corner breakage and particle splitting and b evolution of corner breakage

Full size image

The observed low levels of splitting, both in absolute and relative terms in comparison with corner breakage, may raise the question of whether particle splitting could have been omitted from the breakage model. Traffic loading for \(H\) = 6 m, representative of an intermediate to high embankment showing a substantial level of splitting in comparison with other embankments, was repeated starting from an identical embankment but using a degradation algorithm that did not include particle splitting. This resulted in 2.56 mm settlement, compared with 2.10 mm when including also splitting, as shown in Fig. 10. Settlement for the unbreakable assembly was 2.78 mm. A tentative effect of splitting may therefore be assumed. However, more simulations would be needed to confirm such effect, especially given the levels of dispersion shown in Fig. 10. Furthermore, differences in simulation time were minimal.

The rest of this section is dedicated to provide explanations to singular behaviours presented throughout this section. In this regard, it has been discussed that settlements do not seem to be affected by embankment height (Fig. 10) as could have been anticipated. Figure 15 shows the variation of lateral confinement along the half-track symmetry plane with embankment depth, as is common practice in soil mechanics, for selected embankments. Confinements are obtained by integration of normal contact forces over 1.0 m horizontal strips. At-rest values are also included based on:

$${\sigma }_{\mathrm{h}}\left(z\right)={K}_{\mathrm{o}}\cdot \left[{25+\mathrm{z}\cdot \gamma }_{\mathrm{s}}\cdot (1-{n}_{\mathrm{emb}})\right]+{\Delta \sigma }_{\mathrm{h}}$$
(4)

where \({\sigma }_{\mathrm{v}}\) has been adjusted to account for the permanent loading of 25 kPa, \({n}_{\mathrm{emb}}\) is taken for \(H=10\) m, \({K}_{\mathrm{o}}=0.318\) and \({\Delta \sigma }_{\mathrm{h}}\) can be set roughly to 15 kPa (see below). An average increase of lateral confinement with depth is observed, as expected, translating into higher average values of lateral confinement as was the case in Fig. 8. Confinement on the top two to four layers, i.e. top 1 to 2 m, seems to be rather similar for all the studied embankments. Assuming that rearrangement under traffic concentrates on the top layers (see below), this could explain the apparent lack of influence of average lateral confinement and therefore of embankment height on settlements.

Fig. 15
figure 15

Embankment horizontal lateral stress on half-track symmetry plane after traffic loading for breakable tests. Theoretical at-rest horizontal earth-pressure values after Eq. (4) for \(H=10\) m

Full size image

Regarding theoretical at-rest values, a good agreement with the numerical results is observed for low to intermediate embankments. For higher embankments, numerical results seem to underestimate the at-rest values for the deeper portion of the embankment. This discrepancy for the deeper portion can be explained by the use of friction on vertical walls when applying the permanent load of 25 kPa (see Sect. 3.2). The use of friction can be expected to result in the effect of permanent loading concentrating in the top portion of the embankment, where a good agreement with at-rest values is observed, whereas it should dissipate with embankment depth. Contrary to traffic loading, permanent loading represents a uniformly applied load throughout the whole track rather than exclusively over the studied portion, therefore wall friction should have omitted through the permanent loading phase and posteriorly implemented for traffic loading. In fact, the effect of permanent loading on the at rest values is \({K}_{\mathrm{o}}\cdot 25\) kPa, i.e. roughly 8 kPa. This is approximately the observed average difference between at-rest values and numerical values for the deepest portion of the higher embankments. It should also be observed that for the top layers, an average horizontal value close to 23 kPa is obtained. This is in fact equivalent to the aforementioned 8 kPa due to permanent loading together with roughly 15 kPa representing the overall additional effect of the locked-in stresses \({\Delta \sigma }_{\mathrm{h}}\) due to compaction and previously set to such value in Sect. 4.1. This actually resulted in a good agreement between at rest values and numerical values due to compaction before permanent loading (Fig. 8).

Settlements are originated by particle rearrangement, which can be correlated to particle rotation [10]. The distribution of particle angular velocities at the end of the last unloading cycle is presented in Fig. 16 for selected embankments. Angular velocities are represented by their vector moduli. Distinction between values for the whole embankment and for the top four layers, i.e. approximately top 2 m, is done. Rotational velocities follow rather similar distributions for the top four layers (similarly, results for the top six layers, i.e. top 3 m, still followed near-identical distributions), whereas distinctive variation is observed for the whole embankment values. For these, rotations decrease with embankment height, especially for lower embankments, becoming this decrease minor for higher embankments. This is explained by the fact that the deeper within the embankment, the less rotation is produced by loading of the top layer. Therefore, for high embankments, the more inert part dominates the average values, whereas for low embankment, the top layers still represent a substantial part of the body of the embankment and hence dictate the average values. These results corroborate the above assumption that rearrangement under traffic concentrates on the top layers and hence explains the apparent lack of effect of embankment height on settlements, especially for intermediate to high embankments. It is therefore this substantial scale difference with triaxial conditions that ultimately seems dictate the observed discrepancies on the effect of confinement on settlements.

Fig. 16
figure 16

Distribution of particle rotation velocities for breakable tests: a the whole embankment and b the top four layers

Full size image

5 Concluding remarks

Numerical half-track embankments with heights ranging between 2 and 10 m have been generated by successive dumping and compaction of layers of assemblies of simple tetrahedral clumps with four asperities each. Higher embankments resulted in slightly more compact configurations. Final values of horizontal stresses at the vertical walls seemed to grow with embankment height. A good agreement with theoretical at-rest values was observed, given that the additional effect of locked-in stresses due to compaction was accounted for.

Cyclic loading of the embankments, representing railway traffic, for both breakable and unbreakable assemblies, was carried out. Settlements did not generally seem to follow a discernible pattern with embankment height. However, a clear effect of degradation was observed, largely resulting in increased settlements and in occasions doubling the ones for unbreakable assemblies. This effect was more noticeable in intermediate to high embankments, were uncertainties on settlements development were larger. It is also remarkable that although breakage generally resulted in higher settlements, it also resulted, in occasions, in a more stable logarithmic response with reduced settlements. An analysis of particle rotation velocities and local values of confinement showed that particle rearrangement and hence settlement accumulates mainly in the top layers, thus explaining the above lack of influence of embankment height. Regarding degradation, corner breakage was dominant compared to splitting. Degradation values for embankments were clearly below the ones for triaxial tests from a previous study using the same material model. Furthermore, unlike previous numerical studies for triaxial tests including degradation, given the low levels of breakage, no clear correlation between the evolution patterns of corner breakage and settlement development was found. In fact, settlements did not seem to stabilize or shakedown for the considered number of loading cycles as was the case in previous experimental and numerical studies on railway ballast under triaxial conditions. Instead, logarithmic and exponential growth responses were observed. Regarding the resilient response under cyclic loading, a linear stiffening effect of embankment height was observed with a minor influence of breakage.

The key contribution is thus the modelling of high rockfill embankments including degradation, something not attempted to date, where furthermore the construction process for real embankments is partly accounted for during embankment generation. Among the main conclusions is that the specific scale, boundary and stress conditions of embankments results in a behaviour deviating from that observed for ballast materials under triaxial conditions. Furthermore, the mechanical response is marked by a substantial degree of uncertainty exacerbated by particle degradation. Additionally, common geostatic theories were able to mostly explain the observed results regarding horizontal earth pressures and explanations were provided where deviations were detected.