Dynamic Behavior in Transition Zones and Long-Term Railway Track Performance

Transition zones between embankments and bridges or tunnels are examples of critical assets of the railway infrastructure. These locations often exhibit higher degradations rates, mostly due to the development of differential settlements, which amplify the dynamic train-track interaction, thus further accelerating the development of settlements and deteriorating track components and vehicles. Despite the technical and scientific interest in predicting the long-term behavior of transition zones, few studies have been able to develop a robust approach that could accurately simulate this complex structural response. To address this topic, this work presents a three-dimensional finite element (3D FEM) approach to simulate the long-term behavior of railway tracks at transition zones. The approach considers both plastic deformation of the ballast layer using a high-cycle strain accumulation model and the non-linearity of the dynamic vehicle-track interaction that results from the evolution of the deformed states of the track itself. The results shed some light into the behavior of transition zones and evidence the complex long-term response of this structures and its interdependency with the transient response of the train-track interaction. Aspects that are critical when assessing the performance of these systems are analyzed in detail, which might be of relevance for researchers and practitioners in the design, construction, and maintenance processes.


INTRODUCTION
Railway tracks are known to show higher degradation rates at transition zones between sections with different characteristics and support conditions, i.e., between embankment and bridges, tunnels, or other structures (ERRI, 1999;Banimahd et al., 2012;Varandas et al., 2014a;Paixão et al., 2013;Boler et al., 2018). These locations often exhibit differential settlements that develop during the life cycle of the track due to multiple structural and geo-mechanical reasons. Because the permanent deformation accumulation at transition zone differs significantly from the adjacent locations, differential settlements develop, leading to deteriorated longitudinal level profiles of the track, manifested as a subsidence spanning over a length of a few meters (typically identified in the wavelength range of 3 m < λ < 70 m (CEN, 2017) in the track geometry data obtained by track inspection vehicles). In the long-term, the aggravation of such differential settlements significantly contributes to the deterioration of the track geometric quality and to the amplification of the dynamic vehicle-track interaction, which further increases the deterioration process itself. Eventually, track maintenance interventions are required to reestablish adequate track geometry and safety levels, which is costly, time-consuming, reduces track availability and disturbs normal train operation.
Although many numerical studies have focused on the influence of the vertical stiffness variation to this problem (Banimahd et al., 2012;Costa et al., 2014;Paixão et al., 2018), few have addressed a more critical aspect: the impact and longterm development of differential settlements (Paixão et al., 2016c;Stark and Wilk, 2016). This is probably because it has been far easier to simulate the transient response of a vehicle crossing a transition zone, using adequate numerical tools, than to simulate the long-term behavior of this system, not only because it involves more computational effort, but also because it requires more robust models and mathematical formulations.
Most studies addressing track settlements at transition zones considering the train-track interaction have used simplified 1D Winkler type models coupled with empirical settlement formulas, as in (Mauer, 1995;Hunt, 1997;Kempfert and Hu, 1999;Varandas et al., 2014a). However, because 1D Winkler models require extensive calibration of material parameters and other modeling aspects, such as equivalent springs and dampers (Punetha et al., 2020), its application to real scenarios or to perform parametric studies has been somewhat limited.
While many authors have studied the transient train-track interaction at transition zones using tri-dimensional numerical approaches Paixão et al., 2018;Ngamkhanong et al., 2020), others have proposed 3D FEM models with true tridimensional plastic deformation formulations to simulate the development of settlements, but mostly focusing on aspects other than transition zones. For example, Li et al. (2016) and Shih et al. (2019) presented implementations in commercial FEM software, but limited the analysis to very short sections of track and considered the static loading only (no vehicle-track interaction). Shan et al. (2017) developed a more robust framework to study slab track transition zones to bridges that considered the vehicle-track dynamic interaction to study the impact of the differential settlements at the approach to the structures. Despite the advances, the study only accounted for the plastic deformation of the subgrade, consequently, not applicable to ballasted tracks, where the ballast layer frequently plays a major role in the development of settlements (Selig and Waters, 1994). Moreover, the study adopted an empirical plastic model that was updated at regular load cycle increments, thus, not necessarily considering the amplitude of the plastic strains to update the deformed shape of the track and its impact on the vehicle-track interaction.
Although very robust numerical approaches have been proposed to simulate railway tracks and the complex degradation behavior of the particulate media of the ballast layer, such as the Discrete Element Modeling (Dahal and Mishra, 2020), its application is still limited to very small stretches of track (models with only a few sleepers long) (Chen and McDowell, 2016) because of the extreme computational effort that is required to simulate structures such as transition zones.
To obtain more insight into the long-term behavior of transition zones in ballasted tracks and considering the available tools and the computational challenging task, the authors adopted a full tri-dimensional FEM approach (Varandas et al., 2020) that has the following main features: 1) it incorporates a robust high-cycle strain accumulation model of the ballast layer based on the classical plasticity theory that considers both the frictional and volumetric degradation mechanisms (Suiker and de Borst, 2003); 2) it considers the dynamic vehicle-track interaction using a non-linear hertzian contact formulation, the sleeper-ballast interaction using a bilinear formulation and the consequent non-linearity that results from the vehicles running on uneven longitudinal track profiles; 3) it accounts for the non-linear resilient behavior of the ballast layer, which is critical to adequately evaluate the stress levels that, in turn, are required as inputs for the strain accumulation model.

Dynamic 3D Train-Track System Modeling in FEM
The numerical approach that the authors used in this work resorts to a three-dimensional FEM analysis program to calculate the dynamic train-track interaction, designated Pegasus (Varandas, 2013). It is fully coded in MATLAB ® and has been developed over the past years with the aim to build a FEM code capable to compute, in admissible time with reasonable computational power, the response of a general nonhomogeneous railway track system, to obtain more insight into the behavior of such systems. For this reason, it has often been used to study railway transition zones with irregular track longitudinal profiles (Paixão et al., 2018;Varandas et al., 2017;Varandas et al., 2016b;Paixão et al., 2016d;Varandas et al., 2014b). One of the main advantages of the program is that it can consider both the train-track interaction and the non-linear constitutive behavior of the granular layers supporting the track superstructure Paixão et al., 2016a).
The program establishes three distinct systems that interact via contact forces ( Figure 1): 1) the vehicle model; 2) the track superstructure (rails, sleepers and fastening system); 3) and the ballast-substructure system (the ballast and underlaying soil layers). The vehicle system is represented by rigid masses connected by spring-damper elements. The rails and sleepers of the track superstructure system are represented by Euler-Bernoulli beam elements and the rail pads as spring-damper elements. Given the frequencies under analysis in this work, up to about 80 Hz, and the focus of the paper in settlements predictions, the authors considered that the representation of the rails with the Euler-Bernoulli beam formulation was an adequate approach (Knothe and Grassie, 1993), and that the consideration of a more complex and computationally demanding beam formulation would not yield different longterm results. The ballast-substructure system comprises loworder fully integrated eight-node solid hexahedral elements.
The sleeper-ballast interaction forces in the vertical direction follow a bi-linear formulation to account for the eventual development of unsupported sleepers, or also known as "hanging sleepers" (Varandas et al., 2016b). The wheel-rail interaction forces follow a non-linear hertzian contact formulation given by F k c δ 1.5 , in case of contact, where δ is the indentation of the corresponding wheel and k c is a normal contact stiffness, herein considered constant and equal to 1.0 × 10 8 kN/m 1.5 (Newton and Clark, 1979;Varandas et al., 2017).
Typical analyses with the program comprise a sequence of static and dynamic calculation stages. The Conjugate Gradient Method (Bathe, 1996) is used in the static analyses. A mixed implicit-explicit integration scheme is used in the dynamic analyses to perform the time integration, aiming at reducing the required computational time, as explained in (Varandas et al., 2017): the track superstructure is integrated with the implicit Newmark constant acceleration method, and the ballast-substructure system is integrated with the explicit scheme presented by Zhai (1996).
Regarding the material behavior, linear elasticity is assumed for all materials. An exception is made for the ballast layer which follows a nonlinear elastic behavior for a more accurate simulation. This accounts for the typical change in material stiffness with the stress level, when loaded in elastic or quasielastic regime (Lekarp et al., 2000). This aspect is more relevant in the ballast layer because it undergoes higher stress amplitudes in the cyclic loading than the layers below. The non-linear elastic formulation considered in Pegasus is the pressure dependent hypo-elastic material law commonly known as k − θ model (Brown and Pell 1967). It defines that the resilient modulus, E r , is a function of the sum of the normal stresses, bulk stress, (θ), defined positive for compression, according to Eq. 1: leaving the Poisson's ratio constant. The reference stress, θ 0 , is taken as 100 kPa. The vehicle is modeled as a single bogie, resting on the primary suspension of two wheelsets, and a half carbody mounted on the secondary suspension, as depicted in Figure 2. The vehicle model comprises eight degrees-of-freedom (dofs): four translations and four rotations.

Ballast High-Cycle Strain Accumulation Model
The approach adopted in this study incorporates a high-cycle strain accumulation model to estimate the deformation accumulation process caused by many load cycles imposed by passing vehicles. To avoid calculating the complete (eventually non-elastic) stress/ strain path caused by every single axle pass, which would be impractical and prone to errors (Niemunis et al., 2005), the  process adopted by the authors consists of calculating the envelope of the maximum residual deformations generated over individual load cycles, allowing for the calculation of the total accumulated settlement over a relatively large number of load cycles, following the formulation of the high-cycle strain accumulation model presented by (Suiker and de Borst, 2003). This model was parameterized with triaxial tests on ballast and sub-ballast aggregates (Suiker, 2002) and was developed for the application to railway track simulations. It was selected for this study because it is suited to account for the significant force non-linearity that results from the variable vehicle-track interaction caused by deformed states of the track itself. An aspect worthy of notice is that the accumulation of strain (densification due to cyclic loading) is split in two mechanisms, which may develop simultaneously and independently: 1) volumetric compaction of the material; 2) and frictional sliding of the particles. The approach describes only the accumulation of the average permanent deformation per load cycle. The short-term (dynamic) event is decoupled from the long-term (deformation) event, in which the application of load cycles is implicitly accounted for by the constitutive model. It is described in mechanical stress definition (compression is negative), although the mean normal stress invariant p follows the classical geotechnical definition (i.e., positive in compression). This makes that both invariants, the mean stress, p, and the deviatoric stress, q, are generally positive: The stress due only to gravitational loads is denoted by the subscript "geo", and the stress due solely to the train loading itself, therefore excluding the gravitational component, is denoted by the subscript "cyc": σ cyc (t) σ(t) − σ geo (t). In the presented equations, bold letters are used for tensor/vector quantities. The volumetric and deviatoric strains are respectively defined as: The development of the model is based on the classical plasticity theory (yield surfaces, flow rule, hardening, associated plasticity), in that the yield (or shakedown) surface corresponds to the Drucker-Prager cone, with a compression cap. The permanent strain increment, Δε p , generated in each cyclic loading process is therefore obtained with: where Δε p q is the amplitude of the permanent deviatoric strain increment, Δε p v is the amplitude of the permanent volumetric strain increment and m f and m c are the tensors that define the directions of the deviatoric and volumetric plastic flows, respectively. The superscript p in Eq. 6 refers to permanent, which is omitted for simplicity hereafter because all strains in this accumulation model are permanent.
The amplitude of deviatoric strain increment is obtained by: where α f and c f are model parameters, ( q/p) cyc is the peak q/p ratio, exclusively caused by the vehicle loading itself, and h f (ε q ) is the shakedown evolution function of the frictional mechanism, being ε q the total accumulated deviatoric strain (the history parameter of the frictional mechanism). The amplitude of volumetric strain incremental due to the compaction mechanism only is obtained by: where α c and c c are model parameters, pressure p 0 defines the initial consolidation of the granular material, ( p/p 0 ) cyc is the peak p/p 0 ratio, again considering only the vehicle loading, h c (ε v,c ) is the shakedown evolution function of the compaction mechanism, being ε v,c the total accumulated volumetric strain caused by volumetric compaction (the history parameter of the compaction mechanism). Here are the Macaulay brackets, defined by x x if x > 0 and 0 otherwise. Finally, the incremental amplitude of the volumetric strain is calculated by: In that d f (ε q ) is a function that defines the amount of dilation/ contraction mobilized during deformation by frictional sliding of particles. Following previous calibration works by the research team (Varandas et al., 2020), the parameters of the ballast high-cycle strain accumulation model that were considered for this study are the following: Transition Zone Long-Term Simulations

Description of the Model
A model of a transition zone built with Pegasus is shown in Figure 3, which only depicts the ballast-substructure system,  (Ayachit, 2015). The width of the model is 7.5 m and its length comprises 37 sleepers, spaced 0.60 m, therefore measuring approximately 22.5 m long. The section of the track on embankment is about 12.8 m long and includes a 30 cm ballast layer, a 30 cm sub-ballast layer, resting on top of the subgrade. The section of the track on the bridge deck is about 9.7 m long and only includes the track elements down to the ballast layer. The embankment-bridge interface is located at x 0.3 m. A simplified representation of a wedge-shaped backfill is included in the embankment section, under the ballast layer and adjacent to the bridge end, simulating the design of a typical transition zone with layers in cement bond mixture (Paixão et al., 2015). To properly visualize the mesh of the inner part of the backfill in Figure 3, the elements of the model at x > − 1.75 m and y < 0 m were excluded from the view. The mesh is not uniform along the longitudinal axis, x: the central section (−7 m < x < 4 m) is the region of study and has a finer mesh; the two extremities (x < − 7 m and x > 4 m) were built to ensure stationary moving conditions in the central section and have coarser meshes. The central section comprises 19 sleepers, and the coarser sections comprise nine sleepers each. The model has 197,055 degrees-of-freedom, 65,685 nodes, 1,176 beam elements and 55,832 brick elements.
The geometrical and material properties of the track model are based on data presented in (Paixão et al., 2018), but considering the European standard gauge (1.435 m). The track superstructure is composed of continuously welded UIC60E1 rails, monoblock concrete sleepers (with simplified dimensions 2.6 × 0.3 × 0.19, in m, and weighing 322 kg, each), and Vossloh W14 fastening system with elastomer rail pads Zw700/148/165, with a vertical stiffness value of 160 kN/mm and a damping constant of 9.6 kN/m (Paixão et al., 2018).
The subgrade is assumed to be a homogeneous, isotropic, well compacted, and reasonably stiff soil layer with a depth of 3.7 m.
The upper 0.75 m of the subgrade are represented in the model with brick elements, and the bottom 3 m are represented using an equivalent bi-dimensional viscoelastic Kelvin-Voigt foundation (Varandas, 2013). The lateral walls at the extremities of the model have local dampers to absorb impinging waves (Lysmer and Kuhlemeyer, 1969). The abutment back wall (at x 0.3 m and z < 0.0 m) and the bridge deck surface (at z 0.0 m and x > 0.3 m) were simulated using stiff spring elements in the normal directions to restrain the corresponding displacements. In this study, the bridge and the abutment were considered as rigid bodies, therefore the sagging effect of the bridge deck was not considered. Table 1 summarizes the main properties selected for the track components and the ballast-substructure system. The elastic parameters of the ballast k − θ model were K 1 105 MPa, K 2 0.6 and E min 16 MPa (Aursudkij et al., 2009).
The dynamic simulations in this paper were performed considering the passing of single bogies from the Portuguese tilting passenger vehicle with integrated traction, the Alfa Pendular. The parameters of this bogie model are provided as supplementary material, following the notation presented in Figure 2. The axle load of the Alfa Pendular vehicle is roughly 132 kN. The bogies cross the central region of the model at a speed of 40 m/s (corresponding to 144 km/h), traveling from the embankment to the structure. Considering the track properties, this traveling velocity stays well below the critical speed of the track, and thus the stationary response is quasi-static. Given the high computational effort required to perform the calculations presented herein, the effect of the train speed and the running direction will not be studied in this paper.
Each iteration of the adopted methodology encompasses the gravitational weight calculations, the placing of the vehicle model and the passage of the vehicle. These necessary steps are computationally intensive, mostly because of the consideration of the nonlinear resilient model for the ballast material. Each iteration takes about 50-120 min. to compute, depending on the necessary time step of the integration scheme. It is noted that in every time step of the simulation (in the order of 10 − 5 s) and for every nonlinear element (in this case 14,720 elements) it is also necessary to calculate the sum of the normal stresses θ, the updated secant resilient modulus of the element, E r , using Eq. 1, the updated element stiffness matrix using a fully integration 2 × 2 × 2 Gauss quadrature rule, and corresponding global assembly in the system stiffness matrix.

Long-Term Simulation Methodology
The methodology considered in this study to simulate the longterm behavior of the transition zone consists of a pseudo-time loop approach (Varandas et al., 2014a). The application of such approach to 3D FEM is challenging and a validation of this longterm simulation methodology was presented by the authors in an earlier study (Varandas et al., 2020), focusing on an homogenous foundation scenario. This approach assumes that the effect of the transient response of the train-track system, when trains pass by, can be separated from the long-term plastic deformation behavior of the track bed aggregates and subsequent settlement development, caused by the accumulation of cyclic loads from traffic. Because the plastic deformation of the track bed aggregates is generally several orders of magnitude lower than the resilient deformation, it is an acceptable simplification to perform full elastic train-track dynamic analyses. The methodology also assumes that, in the initial state of the simulation, the track bed layers under the ballast have been heavily compacted during construction and continuously loaded with traffic. As described in detail by Varandas et al. (2020), this iterative approach comprises up to seven steps in each iteration (iter), as depicted in Figure 4. The process starts with a static calculation of an initial state (A) to determine the initial stress field caused by the geostatic load, σ geo . Then, a dynamic calculation is performed (B) to simulate the train-track interaction with the aim to obtain the stress time-history in each finite element caused by the train loading only, σ cyc , followed by the calculation of the stress amplitudes and the peak q p cyc and p p 0 cyc ratios (C), discussed earlier (Eqs. 7, 8).
With that, the corresponding plastic strain increment, Δε p , is calculated (D) with the high-cycle strain accumulation model described earlier (see Ballast high-cycle strain accumulation  (Bowles, 1997;Fortunato et al., 2012;Paixão 2014 model). At this stage, the methodology assumes that successive train passages take place, with unchanged cyclic stress amplitudes, and the plastic strain keeps accumulating, unless a predefined threshold of strain increment is reached, or a specific number of simulated cycles, N, is completed. In the first iteration (iter 1), the methodology also considers the initial compaction conditions (E), which are taken into consideration in the subsequent permanent deformation calculations but do not affect the displacement field. Next, the plastic strain increment Δε p calculated for each element is transferred to the 3D FEM (F) by determining the equivalent element nodal forces (Hughes, 2012). Then, the history parameters of the high-cycle deformation model are updated (G), as well as the shakedown evolution function values of the frictional, h f (ε q ), and compaction mechanisms, h c (ε v,c ). At the end of each iteration (H), with the new deformed configuration of the ballastsubstructure system, a dynamic calculation is performed by placing the track superstructure, which will yield a new track geometry profile. If the total number of simulated cycles is lower than the desired value (I), the process repeats itself from step (B) to calculate the effect of the dynamic train-track interaction with the new track geometry profile, which in turn may alter the dynamic cyclic loading in the ballast-substructure system.

RESULTS AND DISCUSSIONS
In this study the authors present the results of the application of the above methodology to the hypothetical transition zone model presented in Description of the model. It is worth noting that the simulation of this scenario does not correspond to any specific transition zone design. In addition, the authors were only focusing on the impact of the plastic deformation of the ballast layer in the long-term behavior of this system and assumed that the backfill and foundation layers were stabilized. For that reason, the plastic deformation of the subgrade is not accounted for in this study. Although this aspect could be included in this study, it was not the focus of this work. Moreover, this aspect has been addressed in studies by other authors (Shan et al., 2017). In this simulation, the plastic deformation of the ballast layer caused by successive passages of the vehicle Alfa Pendular, from the embankment to the bridge, developed at different rates along the transition zone. This nonhomogeneous behavior caused the growth of differential settlements and the resulting evolution of the longitudinal rail level as depicted in Figure 5. It is visible that the development of the ballast settlements was faster on the structure than on the embankment, where it also stabilized faster -on the embankment, the rail level at N 10E3 is practically identical to that at N 100E3. Due to this, the total settlement in the bridge section is about three times the total settlement attained in the embankment section, although the amplitude of the developed differential settlements at the end of the simulated period is still very limited (less than 0.9 mm) and compatible with normal railway operation. For example, in Europe, the EN 13848-5 (CEN, 2017) establishes that the immediate action limit of isolated defects of the longitudinal level is 23 mm and the recommended alert limit range is 10-17 mm, for speeds between 120 km/h and 160 km/h. In a more realistic scenario, e.g., with even greater number of load cycles, higher amplitudes would eventually been reached, which would require maintenance actions such as tamping and leveling to reestablish adequate track position. Figure 6 shows tri-dimensional views of the deformed shape of the ballast-substructure system obtained at N 100 and at N 100E3, which evidence the concentration of plastic deformation under the sleepers' positions, being more pronounced on the bridge section.
The authors are aware that, in reality, significant settlements may occur in the backfill (caused by the consolidation and densification of the subgrade layers and its foundation) which also contribute to the long-term evolution of the longitudinal rail profile, hence also affecting the dynamic vehicle-track interaction. In fact, the amplitude of the backfill settlements is often greater than that of the scenario simulated here, although this was not the focus of the study presented here, as mentioned earlier. To elucidate on the variability of the evolution of the deformation of the ballast layer, Figure 7 presents examples of settlement accumulation at two positions under two sleepers (z 0.3 m) of the model: 1) at the track center line (CL), y 0 m; 2) under the rail alignment, y 0.75 m. One of the two sleepers is located slightly before the embankment-bridge interface, at x −1.2 m, and the other is located on the structure, at x 0.6 m. The figure presents plots in both linear and log scale with respect to the number of load cycles. The first observation is that the ballast settlement follows a non-linear evolution with the number of load cycles, showing an initial stage of rapid settlement growth, that tends to stabilize after the first 10 thousand load cycles. This is consistent with the shakedown theory (Werkmeister   , 2001). It is also visible in Figure 7 the development of the behavior called "central binding" from an early stage, which is evidenced by the slightly higher amplitude of the ballast settlement under the rail alignment than that on the track center line, leaving the sleeper supported mostly on its central part when the track is at rest, not loaded by the vehicle. This phenomenon is more pronounced under the sleeper on the embankment than under the sleeper on the bridge. Figure 8 presents another aspect of the ballast vertical deformation evolution, assessed in terms of longitudinal settlement profiles of the ballast surface right under the  sleepers (at z 0.3 m), located in the center line of the track (y 0 m). The figure shows that the deformation rates of the sleepers' bedding in the ballast layer away from the embankment-bridge interface is somewhat constant in each side, though at higher rate on the bridge side than on the embankment. However, under the sleepers near the embankment-bridge interface, there are relevant differences in deformation rates between adjacent sleepers that might contribute to the development of hanging sleepers, which will be addressed later. This is a typical problem at transition zones in that some sleepers hang from the rails and their base no longer stays in contact with the ballast layer because the bending stiffness of the track frame (rails, fastening system and sleepers) prevents it from adapting to the uneven surface of the ballast that was created. Irregularities in the rail longitudinal level have strong impact on the dynamic wheel-rail interaction. In fact, this aspect may be considered a key factor to assess the performance of transition zones or other critical locations in railway tracks Mauer, 1995). In view of this, Figure 9 shows the wheelrail interaction forces of both leading and trailing axles, after applying a low pass filter (LPF) with a cut-off frequency of 80 Hz . Apart from a notable force amplitude increase when the vehicle enters the stiffer structure, it is visible, for both axles of the bogie, slight changes in forces with the number of load cycles. For example: 1) there is a slight reduction trend of the forces at the approach to the bridge, between x −2.5 m and x 0 m; 2) and a slight increasing trend right after entering the bridge, between x 0.5 m and x 1.5 m. However, these variations are negligible when compared with other normal and intrinsic causes of variability of the wheel-rail interaction forces, such as the sleeper spacing and the transition from the embankment to the bridge. This negligible effect might be justified by the relatively small settlement that was developed in this simulation, as mentioned above, although it may trigger the development of other problems.
For example, when analyzing the sleeper-ballast interaction forces, it becomes apparent that something more complex is taking place at this level. Figure 10 shows that the resulting sleeper-ballast interaction forces at each sleeper are evolving at different rates along the transition zone: increasing at the approach to the bridge; decreasing around the embankmentbridge interface; and increasing again on the bridge. It is noted that the sleeper-ballast interaction forces were calculated as the sum of all vertical nodal contact forces under each sleeper. The fact that sleeper-ballast interaction forces are decreasing near the interface suggests that hanging sleepers might be developing at that location and that the adjacent sleepers that are not hanging are receiving the additional load that the hanging sleepers are not able to transfer to the ballast layer.
To investigate the development of hanging sleepers, Figure 11 shows the transient displacement undergone by the ballast at different positions, x, along the transition zone, evaluated right under the sleepers, aligned with and below the rail (at y 0.75 m). The red lines (on top) correspond to the first load cycle, while the black lines (at the bottom) correspond to cycle N 100E3. The displacements of the corresponding sleepers are plotted in the same graph using doted lines. It is visible that, in general, the doted lines are overlaid by the full lines, indicating an almost permanent contact between the sleepers' base and the ballast layer. However, due to the differential growth of settlements, some hanging sleepers were developed. This behavior is particularly visible in sleepers at positions x 0 m and x 0.6 m in load cycle N 100E3 because they are not in contact with the ballast before the axles pass by (dotted and full lines do not match) and then they barely touch the ballast when the axles pass directly over it. This behavior is a known malfunction, particularly common in transition zones, that is reported to further promote track degradation (Lundqvist and Dahlberg, 2005;Zhu et al., 2011).
Another aspect that helps to assess the performance of transition zones and the structural consequences of the FIGURE 10 | Evolution of peak vertical sleeper-ballast interaction forces with N.
Frontiers in Built Environment | www.frontiersin.org March 2021 | Volume 7 | Article 658909 development of differential settlements and unsupported sleepers is to analyze the change in track vertical stiffness along the transition zone. Here, the authors calculated the track vertical stiffness as perceived by the axles of the vehicle by calculating the ratio between the wheel-rail dynamic interaction force (presented in Figure 9) and the maximum resilient vertical displacement of the rail due to that axle load (i.e., the amplitudes of the rail displacement depicted in Figure 11). Figure 12 presents the obtained vertical stiffness profiles along the transition zone at different load cycle numbers, where it is visible the typical increase in vertical stiffness as the track approaches the bridge. In this case, the track stiffness increases about 60% when entering the bridge. It is also clear in Figure 12 the relevant change in the profiles with the number of load cycles, especially around the embankment-bridge interface. There is a significant reduction trend in stiffness at this location and it appears that the length in which the stiffness variation takes place increases with N: at N 1, the length of the variation is about 2.5 m and at N 100E3 increases to about 4 m. A slight reduction in stiffness is also apparent on the embankment section, which might be a consequence of the sleepers' central biding phenomenon mentioned earlier.
One of the reasons to develop 3D FEM models to assess the performance of transition zone is that not only it allows for a more accurate reproduction of the load distribution and wave propagation in the structure, but also enables a proper stress analysis inside the ballast and subgrade layers, which is crucial to understand in greater depth the structural and geomechanical phenomena that might be taking place. Herein, the authors will analyze the stress levels in the ballast layer at the relative locations A to D identified in Figure 13, for a given track cross section. Figure 13 also depicts the transient variation of the mean, p, and deviatoric, q, stresses in the ballast layer under the sleepers located at x −1.2 m at x 0.6 m, as the vehicle passed by, at different stages of the cyclic loading. A low-pass filter was also applied to these results with a cut-off frequency of 80 Hz.
The results regarding the elements at positions B and D under the sleeper at x −1.2 m suggest that this position is receiving additional load from adjacent unsupported sleepers because the stress levels at rest are increasing with the FIGURE 11 | Comparison between the deflections of the sleepers (doted lines) and of the ballast layer (full lines) as the vehicles passes by at load cycle N 1 (red) and N 100E3 (black). number of load cycles. However, under the sleeper at x 0.6 m, which is hanging (see Figure 11), the load transfer to the ballast has undergone a reduction with the number of load cycles. Both these observations agree with what was evidenced in Figure 10 about the evolution of the sleeper-ballast interaction forces. Despite the above observation that some locations in the ballast layer seem to be experiencing an increase in the stress level with the evolution in loading cycles, such as under the sleepers located between x −2 m and x 0 m (Figure 10), there are other more relevant factors that govern the development of plastic deformation. In particular, the high-cycle strain accumulation model implemented in this study (see Ballast high-cycle strain accumulation model) considers the ratio ( q/p) cyc as a critical factor for the accumulation of plastic deformation in terms of deviatoric strain, mainly because it indicates a stress level that is closer to the plastic yielding surface. In this formulation, only the cyclic amplitude is accounted for, that is, only the stress increment caused by the vehicle is considered for plastic deformation. Figure 14 presents the obtained peak ratios, assessed at the relative element positions A to D at every sleeper position along the transition zone. A few observations are worth noting. Firstly, it is visible that the elements located lower in the ballast layer (C and D) undergo higher peak q/p ratios although the stress levels, in terms of q and p, are somewhat smaller (Figure 13), compared to what takes place more superficially (i.e., A and B). Secondly, the locations where an increase in sleeper-ballast interaction forces was identified do not seem to correspond to the locations denoting higher peak q/p ratiosthese are clearly located in the track section on the bridge, being the first sleeper on the bridge the worst position at x 0.6 m. Another interesting observation is the similarity between the track stiffness variation profiles presented earlier in Figure 12 and the peak q/p ratios profiles in Figure 14, which suggests a correlation between these two factors. Many studies have been developed trying to establish the cause-effect correlation between track stiffness and track degradation (Pita et al., 2004;Quibel et al., 2010;Powrie and Le Pen, 2016). And lastly, the variation of the amplitude of the peak q/p ratios shown in Figure 14 seems to correlate quite well with the long-term deformed configuration of the system presented earlier (see Figures 5, 8), which is somewhat expected given that the formulation of the model explicitly takes into account the amplitude of that parameter.

FINAL REMARKS
In this work, the authors presented the application of a methodology to simulate the long-term behavior of a hypothetical railway track transition zone scenario aiming at shedding some light into the complex behavior of these structures. In particular, the study focused on the plastic deformation growth of the ballast layer, considering the dynamic train-track interaction, and its interdependent impact on the long-term behavior. The adopted methodology considers a robust high-cycle strain accumulation model to estimate the plastic deformation of the ballast layer and uses a pseudo-time loop approach in combination with a 3D FEM train-tracksubstructure model for dynamic analyses. It was demonstrated the suitability of the tool Pegasus to analyze the evolution of important aspects regarding the performance of transition zones with the number of load cycles, including: the wheel-rail and sleeper-ballast interaction forces; track vertical stiffness; deformation profiles and 3D configurations; transient responses of the system; and in-depth stress analyses. The case study considered herein corresponds to an approach to a bridge with a relatively sharp increase in stiffness (the stiffness increases by about 60% in about 3 m length) and where it is considered that, hypothetically, the subgrade layers in the embankment are fully stabilized, showing no plastic behavior. Only the plastic deformation of the ballast layer is considered. In these conditions, the long-term calculations showed a gradual development of differential settlements in the transition zone, in the sense of greater settlements on the structure than on the embankment, but without an extension that indicates any imminent need for a maintenance intervention. Furthermore, the development of differential settlements in this direction is somewhat contrary to the expected effect in problematic transition zones of railway bridges, which are associated with a lower leveling in the approach sections. This finding points to the conclusion that the transition zones that are typically problematic are mainly due to deficient compaction and stabilization of the soil layers that form the embankment next to the abutment.
The implemented methodology allowed to evaluate the interdependence between dynamic effects, such as the wheelrail interaction, and long-term effects, such as the evolution of the vertical stiffness at the rail level, of the transition zone. The results evidence that the long-term response of these structures is quite complex and is greatly influenced by the transient response of the train-track-substructure interaction, with a clear feedback loop between these two processes.
In particular, it is shown that the development of hanging sleepers is a natural consequence of the expected development of differential settlements and that variations in stress levels in the ballast layer are strongly influenced by changes in the support conditions of the sleepers. The presented results also show that the increase in ballast support stiffness causes a moderate increase in the forces transmitted by the sleepers, and a very significant increase in the q/p ratio, that will accentuate the uneven development of plastic deformation in the ballast layer.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

FUNDING
Part of this work was financially supported by: Base Funding-UIDB/04708/2020 of the CONSTRUCT-Instituto de I&D em Estruturas e Construções-funded by national funds through the FCT/MCTES (PIDDAC).