Hidden rock joints may fool tunnel safety models in the Himalaya, study warns

For engineers planning hydropower tunnels through the steep, earthquake-prone valleys of Nepal’s Lesser Himalaya, the difference between a slope that stands and one that slides may lie in a modeling choice made long before construction begins. A new finite element study of the inlet portal slope of the Sani Bheri 3 Hydroelectric Project in eastern Rukum District shows that when the dominant foliation joints in the local phyllite rock are explicitly included in the model, the predicted safety margin of the slope drops dramatically — in the worst simulated case, from a factor of safety of 2.15 to just 1.41. The work, published open access in Discover Geoscience, argues that simplified modeling approaches commonly used in early design stages can seriously overestimate the stability of structurally complex Himalayan rock slopes.

The research team, led by Nimesh Bhattarai with co-authors including Tulasi Ram Bhattarai, Abhay Kumar Mandal and Keshab Sharma, compared two ways of representing the same slope numerically. The first, an equivalent continuum model, treats the fractured rock mass as a single homogeneous material whose averaged strength parameters implicitly absorb the effects of joints, weathering and fracturing. The second, a combined continuum model, explicitly represents the slope’s most persistent and closely spaced discontinuity — the foliation joint set, oriented at 36 degrees dip toward 196 degrees — using the Barton–Bandis slip criterion, which describes how joint roughness, joint-wall strength and stress-dependent stiffness govern sliding along weak planes. Both models used rock mass properties derived from the generalized Hoek–Brown criterion, with a Geological Strength Index of 45, an intact rock constant of 7 and a disturbance factor of zero, and were run through an identical series of loading scenarios: dry conditions, steady-state groundwater, transient rainfall infiltration, pseudo-static seismic loading, and combinations of these.

The site itself amplifies the stakes. The slope sits in the tectonically active Lesser Himalaya, bounded by the Main Central Thrust to the north and the Main Boundary Thrust to the south, and is composed of foliated phyllite of the Ranimatta Formation — a low-grade metamorphic rock with pronounced anisotropic behavior. Field mapping identified three major joint sets, moderately spaced with persistence generally under five meters, some filled with sand and silt that reduces shear strength and eases water infiltration. Because the foliation is the most persistent and closely spaced of the three sets, it was selected for explicit representation, while the two wider-spaced sets were omitted to keep the two-dimensional model tractable.

Under the baseline dry condition, the gap between the two approaches was already striking. The equivalent continuum model returned a factor of safety of 4.3 with a maximum total displacement of just 0.062 meters, signaling a highly stable, structurally competent slope. The combined continuum model, however, predicted a factor of safety of only 3.34 and a displacement of 0.21 meters — more than three times larger. The authors attribute this to structurally controlled weak planes that facilitate shear slip and localize deformation, reducing the slope’s overall stiffness and strength in ways that homogenized parameters cannot capture.

Adding water widened that gap further. When a groundwater table was introduced at roughly mid-height of the slope, the factor of safety in the equivalent continuum model fell about 13 percent, from 4.3 to 3.73, while in the combined model it fell about 14 percent, from 3.34 to 2.87, with maximum displacement rising to 0.228 meters. Notably, the pore pressure distributions were essentially identical in both models, since the same hydraulic properties and boundary conditions applied, with positive pressures up to 0.51 MPa in the lower slope. The difference lay in where and how that pore pressure acted mechanically. In the equivalent model, elevated pore pressure reduced effective stress uniformly across the rock mass. In the combined model, it acted directly along pre-existing foliation joints where normal stress was already low, cutting effective normal stress across the weak planes and promoting shear slip. Because joint normal stiffness in the Barton–Bandis formulation decreases as effective normal stress declines, the joints tended to open and dilate, producing markedly larger displacements.

Transient rainfall infiltration introduced a time-dependent dimension to the problem. The team simulated two 24-hour rainfall scenarios: a moderate intensity of about 8.2 millimeters per day, corresponding to half the rock’s hydraulic conductivity of 1.9 × 10⁻⁷ meters per second, and a high intensity of about 16.4 millimeters per day, equal to the full conductivity. In the moderate case, the wetting front crept downward over hours, with pore pressure changes confined to near-surface zones during the first six hours before progressively eroding matric suction at depth. In the high-intensity case, the response was much faster: significant pore pressure developed within the first three to six hours, and by twelve hours pore pressure had turned positive throughout the slope, indicating that infiltration was outpacing the rock mass’s drainage capacity. By eighteen to twenty-four hours, near-saturation conditions had stabilized.

Those hydraulic differences translated directly into stability. Under transient infiltration, the equivalent continuum model’s factor of safety dropped roughly 16 percent from its dry baseline, to 3.60, with displacement rising to 0.082 meters. The combined continuum model fell about 22 percent, from 3.34 to 2.62, with maximum deformation reaching 0.24 meters. Two coupled mechanisms drove the degradation: loss of matric suction in the unsaturated zone stripped apparent cohesion from near-surface layers, while rising positive pore pressure at depth slashed effective stress along the discontinuities. Once again, the key distinction was not how much pore pressure built up, but whether it pressed uniformly against a homogeneous mass or pried open genuine planes of weakness.

Seismic loading amplified the divide. Using a pseudo-static approach calibrated to the Nepal Building Code, which specifies a peak ground acceleration of 0.3 g for a 475-year return period, the team applied horizontal and vertical seismic coefficients of 0.3 and 0.15 simultaneously. In dry conditions, the equivalent model’s factor of safety fell about 40 percent, to 2.57, while the combined model’s fell about 44 percent, to 1.86. When steady-state groundwater was added, the equivalent model dropped 47 percent from its dry baseline to a factor of safety of 2.29, while the combined model’s fell 52 percent relative to its own dry condition — a nonlinear escalation of instability that the equivalent approach did not reproduce. The most critical scenario combined transient rainfall with pseudo-static shaking: the equivalent model bottomed out at 2.15, roughly a 50 percent overall reduction, while the combined model reached 1.41, a 58 percent drop, as transient pore pressure and inertial forces interacted along the foliation planes to strip effective normal stress and mobilize shear failure rapidly.

One counterintuitive detail emerged in the equivalent model under combined rainfall and seismic loading: displacement actually decreased relative to pre-seismic conditions. The authors explain this as an artifact of earlier failure onset — when strength is already reduced, the slope reaches its limiting equilibrium at lower deformation levels, so failure occurs before large displacements can accumulate. The lesson for practitioners is that deformation alone can be a misleading indicator of proximity to failure, another reason the researchers caution against relying on simplified screening models for final design decisions.

The findings come with caveats the authors are careful to spell out. The analysis is two-dimensional and plane-strain, so it cannot capture three-dimensional failure modes such as wedge or toppling involving the site’s other two joint sets, which were deliberately excluded after being judged less influential. The groundwater table was assumed rather than measured, the porosity value of 5 percent was adopted from published literature on fractured phyllite because site-specific data were unavailable, and the pseudo-static seismic treatment cannot represent wave propagation, frequency effects, amplification or strain-dependent stiffness degradation, making its results screening-level indicators rather than predictions of actual coseismic behavior. Because the project remains at the pre-construction stage, the numerical predictions could not be validated against field monitoring. The team recommends piezometers, inclinometers and extensometers during construction, along with three-dimensional discontinuum or discrete fracture network modeling, site-specific hydraulic characterization and dynamic time-history seismic analysis for final design verification. Even so, the central message stands: in foliated, tectonically active rock masses, what a stability model leaves out may matter far more than what it puts in.

Subject of Research: Numerical modeling of coupled hydro-seismic rock slope stability at a Himalayan tunnel portal

Article Title: Finite element analysis of coupled hydro seismic effects on rock slope stability at a tunnel portal in the Lesser Himalaya

Article References: Bhattarai, N., Bhattarai, T. R., Mandal, A. K., & Sharma, K. (2026). Finite element analysis of coupled hydro seismic effects on rock slope stability at a tunnel portal in the Lesser Himalaya. Discover Geoscience, 4(1), Article 415. https://doi.org/10.1007/s44288-026-00789-8

Image Credits: AI Generated

DOI: 10.1007/s44288-026-00789-8

Keywords: rock slope stability, finite element analysis, Himalaya, tunnel portal, foliation joints, groundwater, rainfall infiltration, pseudo-static seismic analysis, hydropower, Barton-Bandis, factor of safety, Nepal