Response mechanism of qanat tunnel under high-speed railway subgrade load considering the influence of intersection angle: 3D DEM-FDM numerical study

Abstract

This paper investigates the three-dimensional mechanism of a qanat tunnel when a high-speed railway alignment traverses a qanat region at an oblique angle, a scenario that differs substantially from the two-dimensional (parallel) case and becomes increasingly complex with varying intersection angles. To address this problem, three-dimensional full-scale coupled Discrete Element Method–Finite Difference Method (DEM-FDM) numerical models are established using FLAC3D and PFC3D for four intersection angles between the railway subgrade and the qanat tunnel axis, namely, 0°, 30°, 60°, and 90°, considering both sandy and clayey soils. The ground response is systematically characterized from macro-scale perspectives—including stress distribution, displacement field, tunnel crown settlement, and subgrade centreline settlement—and from micro-scale perspectives, such as inter-particle force chains and particle contact fabric. The key findings are as follows: (1) irrespective of soil type, a larger intersection angle reduces both the extent of subgrade-load influence on the qanat tunnel and the influence of the qanat on subgrade surface settlement; (2) the settlement of the tunnel crown transitions from a uniform distribution at 0° to a characteristic central-peak pattern at larger angles, with the settlements at 60° and 90° being nearly identical; (3) at the micro-scale, the intensity of force chains around the tunnel decreases progressively with increasing intersection angle, while the contact force fabric reveals that the magnitude of vertical contact forces increases with the intersection angle, indicating a higher overall ground bearing capacity; and (4) under identical loading conditions, clayey soil consistently produces smaller displacements and stronger contact force networks than sandy soil, confirming its superior load-resistance capability. These findings provide a three-dimensional mechanistic understanding of qanat–subgrade interaction and offer quantitative reference data for high-speed railway design in qanat regions along the Belt and Road Initiative corridor.

1 Introduction

Qanats are ancient gravity-fed underground water conveyance systems that collect groundwater through a gently sloping tunnel connected to the surface by a series of vertical access shafts (Yazdi and Khaneiki, 2016; Taghavi-Jeloudar et al., 2013). They are widely distributed across the Xinjiang of China, Iran, Afghanistan, and many countries or regains along the Belt and Road Initiative (BRI) corridor and have been listed as key national cultural heritage objects in China (Ebrahimi et al., 2021). With the rapid expansion of high-speed railway (HSR) networks in the BRI region, an increasing number of railway alignments are required to cross qanat zones, as shown in Figure 1. When a railway subgrade is constructed above a qanat, the embankment self-weight and train loads impose additional stress on the ground, potentially destabilizing the qanat tunnel through progressive collapse (Caudron et al., 2006). How to scientifically evaluate and quantify the interaction between high-speed railway embankments and qanats under different intersection angles and further reveal the influence of changes in intersection angles on the response mechanism of qanats, has become a core research focus that urgently needs to be addressed to ensure the safe construction of high-speed railway projects in qanat areas (Zhang et al., 2023a).

FIGURE 1

Previous research on the qanat–subgrade interaction problem has established a foundation of knowledge primarily through two-dimensional (2D) studies. Physical model tests using planar trap-door apparatus and particle image velocimetry (PIV) have revealed the macro–micro evolution of soil arching around qanat tunnels under progressive loading, including the formation, disturbance, and reconstruction stages of the soil arch (Al-Naddaf et al., 2017; Zhang et al., 2025; Zhang et al., 2021). Two-dimensional DEM (Discrete Element Method) simulations using PFC2D have complemented the physical tests by providing particle-scale insight into the force chain redistribution and contact fabric anisotropy changes during tunnel excavation and loading (Abierdi et al., 2020; Liu et al., 2026; Marshall et al., 2009). Stability analyses based on the limit equilibrium theory, and random finite element limit analysis have been used to derive critical burial depths and failure probabilities for qanat tunnels under subgrade loads (Lee et al., 2006; Osman, 2010; Sahoo and Kumar, 2013; Zhang et al., 2023b). Field investigations along the Tehran–Isfahan high-speed railway in Iran have characterized the distribution, geometry, and soil properties of qanats in that corridor (Chen et al., 2022).

However, in practice, railway alignments frequently cross qanat tunnels at oblique angles rather than at exactly 0° (parallel). When the subgrade intersects the qanat tunnel at an angle, the loading condition is inherently three-dimensional and cannot be adequately captured by extending 2D results. The three-dimensional collapse mechanism is not simply a spatial superposition of 2D responses but becomes qualitatively different as the angle varies, because the length of the tunnel exposed to the subgrade footprint, the stress shadow created by the vault, and the soil arching geometry all change with the intersection angle. Studies specifically addressing the effect of the intersection angle on the three-dimensional qanat–subgrade interaction are scarce in the literature, representing a critical gap for the practical design of HSR lines that must cross qanat zones.

To address the three-dimensional effect caused by varying intersection angles, a number of researchers have carried out three-dimensional investigations on buried structures under embankment loading, examining soil arching, load transfer, and deformation around tunnels, pipelines, and culverts (Chen et al., 2026; Lin et al., 2019; Yoo et al., 2025). More recently, Wu et al. conducted a three-dimensional analysis of long-term deformation in heavy-haul railway embankments, offering systematic procedures for material calibration and load application that are relevant to the present work (Wu et al., 2026a; Wu et al., 2026b). Despite these advances, most existing three-dimensional studies have focused on either perpendicular or parallel alignments. The effect of the intersection angle—a critical geometric parameter that governs the spatial extent of load influence, the three-dimensional stress shadow around the tunnel—has received little systematic attention, leaving a significant gap in the understanding of the full three-dimensional response mechanisms.

In response to this deficiency, this paper presents a systematic three-dimensional DEM-FDM coupled numerical investigation of the response mechanism of qanat tunnels under HSR subgrade loads, explicitly considering the effect of intersection angle. The DEM-FDM (Discrete Element Method–Finite Difference Method) coupled framework is selected because it combines the micro-scale particle resolution of DEM, which is essential for revealing soil arching and force chain behaviour near the tunnel, with the computational efficiency of FDM for the subgrade and far-field soil. Full-scale models are established for intersection angles of 0°, 30°, 60°, and 90° in both sandy and clayey soil conditions. The results are analysed from both macro-scale (stress distribution, settlement, displacement field) and micro-scale (force chains, contact fabric) perspectives, yielding a mechanistic picture of how the intersection angle governs the three-dimensional ground response to HSR loading in qanat regions.

2 DEM-FDM coupled methods

2.1 Coupled approach

Full-scale modelling of a high-speed railway subgrade interacting with a qanat tunnel at multiple intersection angles demands a numerical framework that simultaneously resolves (1) the granular, discrete character of the soil near the tunnel where large deformations, soil arching, and force chain effects are dominant, and (2) the continuum behaviour of the subgrade fill over the wider model domain where deformations are relatively small. A purely DEM approach at the engineering scale would require tens of millions of particles, making computation intractable. A purely FDM (continuum) approach cannot capture the inter-particle force redistribution and fabric evolution that govern the qanat response at the microscopic level. The DEM-FDM coupled method overcomes both limitations by assigning the DEM domain to the near-tunnel region and the FDM domain to the subgrade and far-field soil. Figure 2 illustrates the coupled approach in this study. In the coupled framework implemented in FLAC 3D and PFC 3D (Itasca), the DEM and FDM domains are linked at a shared interface. At the beginning of each coupling increment, the interface is identified and the particles potentially in contact with the interface are determined from their positions relative to the interface geometry. The contact force between each particle and the interface wall is computed from the particle–wall overlap using a linear contact model. After each DEM cycle, the updated contact forces are transferred to the FDM boundary nodes through a barycentric interpolation scheme, updating the FDM nodal forces and velocities; the FDM model is then advanced to the next equilibrium state. This iterative exchange ensures force and displacement compatibility across the interface throughout the simulation. The coupling procedure follows the algorithm detailed in references (Xu et al., 2021; Zheng et al., 2022).

FIGURE 2

2.2 Governing equations

In the DEM domain, the motion of each particle is governed by Newton’s second law applied to both translational and rotational degrees of freedom:where m is the particle mass, is the position vector, F_c is the resultant contact force, F_g is the gravitational body force, I is the moment of inertia, ω is the angular velocity, and M_c is the resultant contact moment. Inter-particle contacts are described by a parallel-bond model, which transmits both forces and moments and can represent the cementation behaviour of natural soils through bond cohesion and tensile strength parameters (Gong and Liu, 2017; Zhu et al., 2021).

In the FDM domain, the equations of quasi-static equilibrium are solved incrementally:where σ is the Cauchy stress tensor, ρ is the soil density, g is the gravitational acceleration vector, is the strain rate tensor, and v is the velocity field. The constitutive response in the FDM domain follows the Mohr–Coulomb elastic–perfectly plastic model with parameters calibrated from field and laboratory data (Liu et al., 2026; Yin et al., 2020). The Mohr-Coulomb elastic-perfectly plastic model was selected for the FDM domain because the available laboratory triaxial data define the peak-strength envelope of the representative SC and CL foundation soils, and because this model is widely used in preliminary railway-subgrade design where the main objective is to evaluate the static stress redistribution and deformation pattern. The model assumes non-associated plastic flow and does not include cyclic hardening, stiffness degradation, or pore-pressure accumulation; therefore, the present results should be interpreted as quasi-static response mechanisms rather than long-term dynamic performance predictions.

3 Model setup and micro-parameter calibration

3.1 Soil type and macro-scale strength parameters

The numerical model is based on the in-situ conditions of the Tehran–Isfahan high-speed railway project in Iran, which traverses active qanat zones. Field investigations and laboratory tests identified two representative soil types for the qanat foundation: SC (clayey sand) and CL (lean clay), whose macro-scale strength parameters are summarized in Table 1. In the coupled model, the foundation soils surrounding the qanat tunnel are represented as discrete particles using DEM, whereas the railway subgrade fill is modelled as continuum zones in FLAC3D. The DEM micro-parameters were calibrated to ensure that the numerical assemblies reproduce the target Mohr–Coulomb strength parameters given in Table 1. The FDM subgrade fill was assigned compacted granular-fill properties appropriate for high-speed railway embankment analysis; these subgrade parameters are also listed in Table 1.

3.2 DEM micro-parameter calibration via triaxial simulation

Before constructing the coupled model, the DEM micro-parameters for both soil types were calibrated by simulating conventional triaxial compression tests. The triaxial numerical model uses PFC3D particles enclosed within rigid top and bottom loading platens (modelled as walls) and a flexible rubber membrane (modelled by FLAC shell elements in the FLAC-PFC coupled mode), faithfully reproducing the laboratory test boundary conditions (Figure 3). Four confining pressures were applied: 100, 200, 300, and 400 kPa, as shown in Figure 4). The micro-parameters were iteratively adjusted until the simulated peak deviatoric stress at each confining pressure matched the laboratory values with coefficients of determination R² = 0.9998 (sandy soil) and R² = 0.9796 (clayey soil). The linear regression fits confirm a Mohr–Coulomb failure envelope of the form q = 0.8012σ₃ + 5 for sandy soil (φ = 40°, c = 5 kPa) and q = 0.3627σ₃ + 30 for clayey soil (φ = 22°, c = 30 kPa), as shown in Figure 5. The calibrated micro-parameters are summarized in Table 2. The parallel-bond contact model was adopted because both representative foundation soils show measurable apparent cohesion. For the SC soil, the bond strength is intentionally low and serves as a numerical representation of the weak cementation and clay fraction in clayey sand rather than a literal grain-scale cement bond. The DEM particles in the full-scale model are therefore numerical representative particles, not real soil grains. Coarse graining is controlled by macro-scale calibration: the particle assembly is accepted only when the simulated triaxial peak strengths reproduce the target Mohr-Coulomb envelopes. Consequently, the force-chain results are interpreted comparatively among cases with identical numerical resolution and boundary conditions, rather than as absolute real-grain contact forces.

FIGURE 3

FIGURE 4

FIGURE 5

3.3 Model geometry and simulation cases

The coupled model represents a full-scale field cross-section with a model domain of 90 m (length) × 50 m (width) × 8 m (height) for the DEM ground zone, and a FDM subgrade zone of 50 m × 21.6 m × 4 m subdivided into 30 × 10 × 20 structured zones using the Mohr–Coulomb constitutive model. The DEM domain contains approximately 500,000 particles generated by the rain-fall method at a relative compaction of 0.8 (porosity = 0.375). Four intersection angles between the railway centreline and the qanat tunnel axis are modelled: 0° (subgrade parallel to tunnel), 30°, 60°, and 90° (subgrade perpendicular to tunnel). The qanat tunnel has a circular cross-section with an internal diameter D = 2 m and a burial depth of H = 4 m (measured from the ground surface to the tunnel crown), giving a cover-to-diameter ratio H/D = 2. The full simulation cases are summarized in Table 3. The global coordinate system is defined as follows: X is parallel to the railway subgrade centreline, Y is horizontal and perpendicular to the subgrade centreline, and Z is vertical, positive upward. The local coordinate Q is measured along the qanat tunnel axis from the crossing midpoint. For an intersection angle theta between the railway centreline and the qanat axis, Q = X cos(θ) - Y sin(θ). This convention is used consistently for the slicing positions and settlement profiles.

3.4 Model generation procedure

FIGURE 6

3.5 Model validation

The numerical reliability was checked at three levels. First, the DEM micro-parameters were calibrated using triaxial simulations at confining pressures of 100, 200, 300, and 400 kPa, producing coefficients of determination of R² = 0.9998 for SC and R² = 0.9796 for CL. Second, the DEM-FDM coupling procedure was benchmarked against published coupled PFC 3D-FLAC 3D embankment simulations, and the reproduced settlement profile differed by less than 8% at the monitoring points. Third, the 0° parallel-crossing case was compared with the companion two-dimensional physical-model observations for the same cover-to-diameter ratio; The predicted crown settlement of approximately 36 mm is close to the reported experimental value of approximately 34 mm. These comparisons support the use of the model for comparative mechanism analysis, although direct field validation remains a necessary future step.

4 Results

4.1 Effect of intersection angle on stress distribution

Figure 7 presents the plan-view vertical earth pressure distribution at the bottom of the DEM ground domain (z = 0) for the four intersection angles in sandy soil under combined subgrade and train loading. The color scale is consistent across all panels, enabling direct comparison. In all four cases, the earth pressure is elevated within the subgrade footprint area due to the embankment and train loads. The qanat tunnel produces a localized zone of reduced earth pressure along its axis, because the open void disrupts the vertical stress transmission and redistributes stress laterally toward the tunnel walls. The extent and intensity of this stress shadow depend systematically on the intersection angle: At 0°, the tunnel lies entirely beneath the centreline of the subgrade. The stress shadow extends along the full projected length of the tunnel within the subgrade influence zone, and the earth pressure in the centre of the model is notably reduced and relatively uniform. At 30°, the earth pressure above the portion of the tunnel within the subgrade footprint is higher than above the portions outside the footprint, creating an asymmetric pattern. The stress shadow rotates obliquely with the tunnel axis. At 60°, the influence zone of the subgrade on the tunnel is further reduced because only a shorter projected length of the tunnel intersects the subgrade footprint. The stress shadow is noticeably smaller and more localized. At 90°, the subgrade crosses the tunnel at a right angle, and the tunnel influence zone on the stress field is the smallest of all four cases, confined to a narrow band near the tunnel mid-span. These results demonstrate that a larger intersection angle is beneficial for protecting the qanat tunnel, as it reduces the length of the tunnel directly exposed to the elevated subgrade stress, and thereby decreases the risk of tunnel crown settlement and collapse.

FIGURE 7

4.2 Effect of intersection angle on deformation field

Figure 8 presents the plan-view vertical settlement at the bottom of the ground domain for the four intersection angles in sandy soil under full loading. The maximum settlement at the 0° angle reaches approximately 37 mm, occurring in the central region of the model directly above the tunnel centreline. As the intersection angle increases: At 0°, the high-settlement zone is distributed symmetrically along the tunnel axis beneath the subgrade, with the peak occurring at the model centre. At 30°, the region of maximum settlement is oriented at 30° to the Y-direction, reflecting the oblique alignment of the tunnel relative to the subgrade; the settlement peak remains at the model centre (directly beneath the subgrade centreline) and decreases slightly compared to the 0° case. At 60° and 90°, the high-settlement zone contracts further toward the tunnel mid-span, and the maximum settlement magnitudes for these two angles are nearly identical, differing by less than 1 mm. This finding implies that the intersection angle effect on the ground settlement distribution saturates at approximately 60°: increasing the angle from 60° to 90° produces a negligible additional reduction in peak settlement. In all settlement figures and discussions, positive settlement denotes downward movement (subsidence). Therefore, larger positive values indicate greater downward displacement of the subgrade or tunnel crown.

FIGURE 8

Figure 9 shows the vertical settlement of the subgrade surface along the centreline (X-direction) for the four intersection angles. At 0°, the subgrade centreline settlement is largest, reaching approximately 72 mm. As the intersection angle increases to 30°, 60°, and 90°, the centreline settlement decreases monotonically. The settlement profiles are all convex (maximum at the centre, decreasing toward the edges), and the profiles at 60° and 90° are nearly identical, again confirming the saturation of the intersection angle effect beyond 60°. The decrease in subgrade surface settlement with increasing intersection angle is attributable to two coupled effects: (1) a larger intersection angle reduces the exposed length of the tunnel beneath the subgrade, thereby reducing the volume of the stress shadow and the associated differential settlement at the surface; and (2) the more localized stress concentration near the tunnel mid-span at high angles limits the settlement influence zone at the subgrade surface.

FIGURE 9

Figure 10 presents the vertical settlement of the qanat tunnel crown plotted along the tunnel axis (X-direction), with X = 0 corresponding to the mid-span of the subgrade crossing. At the 0° intersection angle, the crown settlement profile is approximately uniform along the full exposed tunnel length, with a peak value of approximately 36 mm, reflecting the fact that the entire tunnel beneath the subgrade is equally affected by the load. As the intersection angle increases to 30°, 60°, and 90°, the crown settlement profile transitions from a flat distribution to a characteristic convex shape (larger at the centre, smaller at the ends), with progressively smaller peak values. The peak values at 60° and 90° are nearly the same and are markedly smaller than those at 0° and 30°. This transition from a flat to a convex settlement profile with increasing intersection angle reflects the progressive reduction of the exposed tunnel length: at high angles, only the section of the tunnel nearest the subgrade centreline is significantly affected by the load, while the tunnel sections farther from the centreline are well outside the subgrade stress influence zone.

FIGURE 10

To characterize the internal deformation pattern at the model centre and to identify the shape of the iso-displacement contours, cross-sections perpendicular to the subgrade direction are extracted at the model midpoint for the four intersection angles (Figure 11). The displacement magnitude is plotted as a function of depth and transverse position. At all four intersection angles, the maximum displacement occurs at the model centre at depth level and diminishes with increasing depth, consistent with the downward propagation of the subgrade load. The iso-displacement contours exhibit a characteristic “V”-shaped pattern (wider at the surface and narrower at depth) at 0° and 30°, where the tunnel void creates a concentrated deformation zone above it. At 60°, the V-shape becomes more rounded (arc-shaped) at the apex. At 90°, the V-shape transforms into a U-shape, indicating that the tunnel influence on the displacement field is more diffuse and symmetric at the perpendicular angle. This evolution of the iso-displacement shape with angle reflects the changing three-dimensional geometry of the stress redistribution zone around the tunnel crown.

FIGURE 11

4.3 Micro-scale analysis

To reduce reliance on visual interpretation, the force-chain discussion is supplemented by quantitative descriptors: the mean normal contact force, the strong-force-chain ratio Rsc (contacts with normal force greater than 1.5 times the mean normal force divided by the total number of contacts), the second-order contact fabric tensor, the fabric anisotropy coefficient, and the vertical-to-horizontal contact-force ratio. These indices are used as comparative measures among simulations with the same particle resolution and post-processing threshold. To investigate the three-dimensional distribution of inter-particle contact forces, cross-sections perpendicular to the tunnel axis are extracted at multiple positions along the tunnel: at Y = 0, 5, 10, 15, 20, and 24 m for the 0°, 30°, and 60° cases, and at X = 0, 5, 10, 15, 20, and 24 m for the 90° case (Y = 0 and X = 0 correspond to the subgrade centreline), as shown in Figure 12.

FIGURE 12

Figure 13 shows the resulting force chain images. At 0°, the force chain distribution is nearly identical across all cross-sections, with a dominant trapezoidal strong-chain zone (indicated by a red frame) in the central region of each section, reflecting the fact that the entire tunnel beneath the subgrade is uniformly loaded. The strong chains are oriented primarily vertically, consistent with the downward propagation of the embankment and train loads. At 30°, the Y = 0 cross-section shows a symmetric strong-chain zone centred on the tunnel, indicating maximum loading influence at this position. As Y increases (moving away from the subgrade centreline), the strong-chain zone progressively migrates away from the tunnel axis and weakens, demonstrating a gradual reduction in the tunnel’s load exposure. The force chains also rotate: at Y = 0 they are predominantly vertical; at Y = 24 m they become predominantly horizontal, reflecting the dominance of the self-weight arch rather than the externally applied load. At 60°, the force chains near the tunnel are noticeably sparser than at 0° and 30°, indicating a smaller loading influence. At 90°, the influence of the subgrade load on the tunnel is most localized: beyond X = 15 m, the force chain cross-sections are essentially identical to the pre-loading state, confirming that the tunnel outside the subgrade centreline zone is unaffected by the train loading.

FIGURE 13

The 30° case is examined in detail by comparing the particle displacement field and the force chain pattern at the same cross-sections (Y = 0, 5, 10, 15, 20, 24 m). Figure 14 presents the displacement (left) and force chain (right) images. From the displacement images: at Y = 0, the displacement field is symmetric about the tunnel centreline, with the maximum displacement at the model centre and a characteristic V-shaped iso-displacement pattern indicating that the tunnel acts as a deformation sink. As Y increases, the maximum displacement zone migrates laterally away from directly above the tunnel, and the tunnel crown displacement decreases progressively. Beyond Y = 10 m, the displacement field transitions from a V-shaped pattern focused on the tunnel to an arc-shaped pattern expanding from the surface, indicating that the load influence on the tunnel is weakening. At Y = 24 m, the tunnel crown displacement is very small (less than 10 mm), showing that the subgrade loading influence extends approximately 20 m laterally from the centreline for the 30° case. From the force chain images: at Y = 0, the force chains around the tunnel are predominantly strong chains, reflecting the high contact forces under full subgrade loading. As Y increases, the proportion of weak chains above the tunnel crown increases progressively, and the force chain orientation rotates from primarily vertical (at Y = 0) to increasingly horizontal (at Y = 24 m). This directional evolution of the force chains indicates a progressive stress redistribution: at positions far from the subgrade centreline, the tunnel is principally under the influence of the self-weight arch (horizontal arch mechanism) rather than the external vertical load, providing a naturally favorable stabilizing condition.

FIGURE 14

To quantify the loading influence zone along the tunnel axis, force chains are extracted at multiple cross-sections perpendicular to the tunnel axis at positions Q = 0, 6, 10, 14, 16, 20 m (where Q is measured along the tunnel from its midpoint). Figure 15 compares the force chain cross-sections for all four intersection angles. The symbol Q represents the position along the qanat tunnel axis. At 0°, the force chain cross-sections are nearly identical at all Q positions, confirming that the loading influence is uniformly distributed along the full tunnel length. In all cases, the force chains on the two sides of the tunnel (lateral walls) are noticeably thicker than those above and below the tunnel (crown and invert), indicating that the lateral walls provide the primary supporting resistance to the load-induced tunnel deformation, while the crown and invert receive comparatively lower contact forces. At 30°, force chains in the range −14 m < Q < 14 m are noticeably stronger and denser than at Q ≥ 16 m. Beyond Q = 16 m, the force chains decrease markedly, indicating that the portion of the tunnel outside the range −14 m < Q < 14 m is largely unaffected by the subgrade load. At 60° and 90°, the strong-chain influence zone is reduced to approximately −10 m < Q < 10 m. These results yield quantitative definitions of the subgrade loading influence length for each intersection angle, summarized in Table 4.

FIGURE 15

4.4 Comparison between sandy and clayey soil conditions

Figure 16 compares the tunnel crown settlement profiles for both soil types at all four intersection angles. The following consistent trends are observed: For both soil types, the crown settlement is largest at 0° and decreases progressively with increasing intersection angle. At 60° and 90°, the crown settlement profiles are nearly identical, confirming the saturation effect of the intersection angle. The qualitative shape of the settlement profile transitions from flat (0°) to convex (30°, 60°, 90°) for both soil types. The quantitative difference between the two soil types is most pronounced at the tunnel mid-span (X = 0) and diminishes toward the tunnel ends (X = ±30 m). At 0°, the sandy soil crown settlement is approximately 4 mm larger than the clayey soil at the mid-span. At higher intersection angles, the differential is smaller at the tunnel mid-span but negligible at the tunnel ends. This differential behaviour arises because the mid-span of the tunnel is most directly affected by the applied load, and this is where the difference in soil cohesion has the greatest effect on the load-carrying ability of the tunnel crown.

FIGURE 16

Figure 17 compares the subgrade centerline settlement on sand and clay foundation. The distribution pattern of subgrade centerline settlement for cohesive soil is similar to that of sandy soil. The displacement of the subgrade is the largest at a 0° intersection angle, approximately 58 mm. As the intersection angle increases, the settlement gradually decreases, with the maximum settlement occurring at the center and less on both sides. At a 0° intersection angle, the subgrade centerline settlement of cohesive soil is smaller than that of sandy soil, by approximately 14 mm. The settlements of cohesive soil at 30°, 60°, and 90° intersection angles are also smaller, indicating that a cohesive soil foundation is more conducive to reducing subgrade settlement.

FIGURE 17

At the 30° intersection angle, the displacement field cross-sections and force chain images for sandy and clayey soil at the same slice positions (Y = 0, 5, 10, 15, 20, 24 m) are compared in Figure 18. The displacement evolution with Y-position is qualitatively identical between the two soil types: as the cross-section moves away from the subgrade centreline, the tunnel crown displacement decreases, and the load influence progressively weakens. The key quantitative difference is that the displacement magnitudes are consistently smaller in the clayey soil cross-sections, confirming higher load-resistance capability.

FIGURE 18

Figure 19 compares the force chains around the qanat tunnel at the Y = 0 section, taking a 30° intersection angle as an example. The raw contact counts are 62,170 for CL and 61,662 for SC; the difference is only approximately 0.82%, so contact number alone is not used as evidence of a stronger network. Instead, the interpretation is based on the combined force-chain pattern, contact-force intensity, and normalized force-chain indices. Under identical loading and numerical resolution, the CL case shows smaller displacement and a more continuous load-transfer path along the tunnel sidewalls, indicating a stronger resistance to deformation than the SC case.

FIGURE 19

5 Discussion

5.1 Physical interpretation of the intersection angle effect

The consistent reduction in ground response severity with increasing intersection angle across all metrics (stress shadow extent, tunnel crown settlement, subgrade surface settlement, force chain intensity, and load influence zone length) can be attributed to a fundamental geometric effect: as the intersection angle increases from 0° to 90°, the projected length of the qanat tunnel beneath the subgrade footprint decreases from a maximum (full tunnel length) at 0° to a minimum (minimum projected length in the limiting perpendicular crossing geometry) at 90°. This reduces the effective span of the void exposed to the subgrade load, thereby limiting the area of disturbed soil arching and the volume of the stress shadow. The consistent reduction in ground response severity with increasing intersection angle across all metrics (stress shadow extent, tunnel crown settlement, subgrade surface settlement, force chain intensity, and load influence zone length) can be attributed to a fundamental geometric effect: as the intersection angle increases from 0° to 90°, the projected length of the qanat tunnel beneath the subgrade footprint decreases from a maximum (full tunnel length) at 0° to a minimum (zero projected length in the extreme theoretical case) at 90°. This reduces the effective span of the void exposed to the subgrade load, thereby limiting the area of disturbed soil arching and the volume of the stress shadow.

The near-identical behaviour at 60° and 90°—where settlements and force chain distributions converge to within 1 mm—suggests that the geometric effect saturates at approximately 60°. Beyond this angle, the projected tunnel length is short enough that the tunnel void is essentially a point load rather than a line load relative to the subgrade width, and further angular rotation does not significantly change the load exposure geometry. This saturation behaviour is practically important: it implies that, when designing HSR alignments in qanat regions, achieving an intersection angle of approximately 60° or greater provides near-optimal protection for the qanat without requiring the more constrained 90° crossing geometry.

5.2 Implications for HSR design in qanat regions

5.3 Limitations and future work

The present model uses an equivalent quasi-static train load and therefore does not capture cyclic stress accumulation, train-speed-dependent dynamic amplification, or long-term degradation under repeated axle passages. In addition, the DEM particles are coarse-grained numerical representative particles. Future work should incorporate moving dynamic train loads, cyclic constitutive behaviour, and field monitoring data from qanat-crossing railway sections to further validate and extend the design implications.

6 Conclusion

References

• 1

  Abierdi, XiangY.ZhongH.GuX.LiuH.ZhangW. (2020). Laboratory model tests and DEM simulations of unloading-induced tunnel failure mechanism. Comput. Mater. Continua63 (2), 825–844. 10.32604/cmc.2020.07946

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  Al-NaddafM.HanJ.JawadS.AbdulrasoolG.XuC. (2017). Investigation of stability of soil arching under surface loading using trapdoor model tests. ICSMGE 2017 - 19th International Conference on Soil Mechanics and Geotechnical Engineering, 2017-Septe (June 2019), 889–892.

  • Google Scholar
  • View reference in article

• 3

  CaudronM.EmeriaultF.KastnerR.MarwanA. (2006). “Collapses of underground cavities and soil-structure interactions _ experimental and numerical models,” in Proceedings of the Proceedings of the 1st Euro Mediterranean Symposium on Advances on Geomaterials and Structures, 3–5.

  • Google Scholar
  • View reference in article

• 4

  ChenW.MinS.ChalaA. T.ZhangY.LiuX. (2022). Assessing compaction of existing railway subgrades using dynamic cone penetration testing. Proceedings of the Institution of Civil Engineers: Geotechnical Engineering175 (4), 439–450. 10.1680/jgeen.19.00303

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  ChenC.YangY.McDowellG.ZhongJ.WangJ.WangL. (2026). DEM–FDM coupled analysis of composite element test for geogrid-reinforced ballast under cyclic loading. Railw. Eng. Sci.2026, 1–17. 10.1007/S40534-026-00438-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  EbrahimiA.MehrabanY.OmidvarbornaH.VakilinejadA.Al-SayighA. R. S. (2021). Kariz (Ancient Aqueduct) system: a review on geoengineering and environmental studies. Environ. Earth Sci.80 (6), 1–13. 10.1007/s12665-021-09545-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  GongJ.LiuJ. (2017). Mechanical transitional behavior of binary mixtures via DEM: effect of differences in contact-type friction coefficients. Comput. Geotechnics85, 1–14. 10.1016/j.compgeo.2016.12.009

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  LeeC. J.WuB. R.ChenH. T.ChiangK. H. (2006). Tunnel stability and arching effects during tunneling in soft clayey soil. Tunn. Undergr. Space Technol.21 (2), 119–132. 10.1016/j.tust.2005.06.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  LinX. T.ChenR. P.WuH. N.ChengH. Z. (2019). Three-dimensional stress-transfer mechanism and soil arching evolution induced by shield tunneling in sandy ground. Tunn. Undergr. Space Technol.93 (March), 103104. 10.1016/j.tust.2019.103104

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  LiuY.LeiH.HuangH.WangM.ManJ. (2026). Face instability mechanisms of shield tunnel undercrossing an existing tunnel: insights from centrifuge model tests and FDM-DEM simulations. Transp. Geotech.60, 102016. 10.1016/j.trgeo.2026.102016

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  MarshallA. M.ElkayamI.KlarA. (2009). “Ground behaviour above tunnels in sand - DEM simulations versus centrifuge test results,” in Euro:Tun 2009 - 2nd International Conference on Computational Methods in Tunnelling, 183–190.

  • Google Scholar
  • View reference in article

• 12

  OsmanA. S. (2010). Stability of unlined twin tunnels in undrained clay. Tunn. Undergr. Space Technol.25 (3), 290–296. 10.1016/j.tust.2010.01.004

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  SahooJ. P.KumarJ. (2013). Stability of long unsupported twin circular tunnels in soils. Tunn. Undergr. Space Technol.38, 326–335. 10.1016/j.tust.2013.07.005

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  Taghavi-JeloudarM.HanM.DavoudiM.KimM. (2013). Review of ancient wisdom of Qanat, and suggestions for future water management. Environ. Eng. Res.18 (2), 57–63. 10.4491/eer.2013.18.2.057

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  WuW.ChenY. T.WanatowskiD. (2026a). Long-term deformation mechanisms of double-track heavy-haul railway embankment under bidirectional shear loadings. Acta Geotech.2026, 1–16. 10.1007/S11440-026-03064-9

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  WuW.ChenY. T.WanatowskiD.WangJ. (2026b). Experimental investigation of heavy haul railway embankment deformation under traffic-induced bi-directional principal stress rotation. Transp. Geotech.60, 102031. 10.1016/J.TRGEO.2026.102031

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  XuZ.ZhangL.PengB.ZhouS. (2021). DEM-FDM numerical investigation on load transfer mechanism of GESC-supported embankment. Comput. Geotechnics138 (June), 104321. 10.1016/j.compgeo.2021.104321

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  YazdiA. A. S.KhaneikiM. L. (2016). Qanat Knowledge: Construction and Maintenance. Springer.

  • Google Scholar
  • View reference in article

• 19

  YinZ. Y.WangP.ZhangF. (2020). Effect of particle shape on the progressive failure of shield tunnel face in granular soils by coupled FDM-DEM method. Tunn. Undergr. Space Technol.100 (April), 103394. 10.1016/j.tust.2020.103394

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  YooM.KwakC.LeeH. U.ParkJ. (2025). Development of seismic fragility curve for subsea railway tunnel based on 3D numerical analysis. KSCE J. Civ. Eng.29 (7), 100149. 10.1016/J.KSCEJ.2024.100149

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  ZhangZ.TaoF.-J.HanJ.YeG.-B.ChengB.-N.XuC. (2021). Arching development in transparent soil during multiple trapdoor movement and surface footing loading. Int. J. Geomechanics21 (3), 04020262. 10.1061/(asce)gm.1943-5622.0001908

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  ZhangY.LiuX.YuanS.SongJ.ChenW.DiasD. (2023a). A two-dimensional experimental study of active progressive failure of deeply buried Qanat tunnels in sandy ground. Soils Found.63 (3), 101323. 10.1016/j.sandf.2023.101323

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  ZhangY.LiuX.YuanS.ZhangT.SongJ.ChenW. (2023b). Probabilistic stability analysis of qanat tunnels in c-φ soil considering soil spatial variability. Eur. J. Environ. Civ. Eng.27 (12), 3763–3783. 10.1080/19648189.2022.2152101

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  ZhangC.LiuX.YuanS.ZhangY.ChenW.WangS. (2025). Critical Depth of Qanat Tunnel Considering Its Interaction with Railway Subgrade: Experimental and Numerical Study. Geomechanics and Geoengineering. 10.1080/17486025.2025.2580001

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  ZhengJ.PrevitaliM.KnappettJ.CiantiaM. O. (2022). “Coupled DEM-FDM investigation of centrifuge acceleration on the response of shallow foundations in soft rocks,” in 10th International Conference on Physical Modelling in Geotechnics. Seoul: Korean Geotechnical Society, 264–267.

  • Google Scholar
  • View reference in article

• 26

  ZhuY.GongJ.NieZ. (2021). Shear behaviours of cohesionless mixed soils using the DEM: the influence of coarse particle shape. Particuology55, 151–165. 10.1016/j.partic.2020.07.002

  • CrossRef
  • Google Scholar
  • View reference in article