Rainfall-induced instability of mountainous photovoltaic slopes under spatially non-uniform infiltration

Rainfall-induced slope failures are among the most frequent and destructive geohazards in mountainous regions. In photovoltaic (PV) installations, panel shading and runoff concentration cause highly non-uniform rainfall infiltration, which alters pore-pressure distribution and threatens slope stability. This study develops a coupled hydro-mechanical modeling framework to investigate slope instability under spatially heterogeneous rainfall. By integrating the Richards equation, the Van Genuchten model, and the strength reduction method within a customized FLAC3D platform, the transient evolution of pore-water pressure, matric suction, and shear strength was simulated for both uniform and non-uniform rainfall conditions. The results show that PV panel shading results in pronounced hydrological inhomogeneity, creating asynchronous suction dissipation and inhomogeneous pore pressure build-up across the slope. Non-uniform infiltration accelerates localized saturation beneath exposed infiltration zones while delaying wetting under shaded areas, resulting in lateral hydraulic gradients and early pore-pressure concentration near the slope toe. The coupling between non-uniform infiltration and interlayer hydraulic contrast is the dominant process of slope failure induced by the PV system. The simulation study illustrates the results caused by such neglect of rainfall heterogeneity can lead to overly conservative estimation of slope stability and inadequate predictions of the localized risk. For mountain photovoltaic construction, the proposed framework can serve as a reliable risk-assessment tool and early-warning design tool in nature equipped with spatially variable rainfall.

1 Introduction

Rainfall-induced slope failures are one of the most widespread and devastating types of geological hazard throughout the world (Dai and Lee, 2002; Crosta and Frattini, 2008; Zhou and Sun, 2019). They exhibit great frequency and pose serious risk to infrastructure and human lives in mountainous areas. More studies have shown that intense rainfall infiltration can affect pore-water pressure and matric suction spatially, which may eventually lead to the reduction in strength, redistribution of stress and shallow landslides, deep landslides (Li et al., 2025; Cao et al., 2025; Li et al., 2026; Li et al., 2024; Zhang SQ. et al., 2025; Cho and Lee, 2001; Iverson, 2000; Fredlund D. G. and Rahardjo H., 1993; Nie et al., 2023). China’s southwestern and southern hilly areas (with steep terrain and most concentrated rainfall) experience rainfall-induced landslide disasters, which account for more than 70% of all geologic disasters recorded in China (Zhang et al., 2014). To understand rainfall-induced instability, many studies have investigated the coupling mechanism between rainfall infiltration and slope stability from different perspectives. Early observation of the field and laboratory (Tu et al., 2009; Wu et al., 2015a; Yuan et al., 2023) phenomena suggested that soil suction rapidly decreases and the pore-water pressure quickly increases, and an active zone was formed at about 1–3 m below the soil surface. These migration and lags of soil moisture will control the onset and evolution of a possible sliding plane. Further research results from Rahimi et al. (2010) and Oh and Lu (2015) show that the dissolution of soil matric suction during rain is one of the crucial causes of slope failure.

Meanwhile, numerical and theoretical studies based on the Van Genuchten model and the Richards equation (Lu and Likos, 2004; Borja and White, 2010) have provided a sound framework for describing the unsaturated-to-saturated transition and corresponding seepage–strength coupling. Building upon these foundations, recent research has introduced more sophisticated hydro-mechanical models to capture the coupled effects of rainfall infiltration, soil structure, and deformation under transient conditions. Subedi et al. (2025) developed a rainfall–seepage–instability coupling approach that links infiltration dynamics with time-varying rainfall, enhancing the understanding of slope response under natural rainfall sequences. Garcia Aristizaba and Riveros Jerez, 2011 further revealed that rainfall redistribution can alter the microstructural pore network of soil, amplify localized infiltration and accelerate suction loss. Similarly, Fan et al. (2015) and Zhou et al. (2024) emphasized the role of antecedent moisture and rainfall sequencing in controlling failure thresholds, while Wang et al. (2025) demonstrated that integrating hydrological modeling with three-dimensional strength reduction analysis significantly improves the predictive accuracy for rainfall-triggered slope failures. Collectively, these findings highlight the importance of characterizing the spatiotemporal variability of rainfall and soil hydraulic response within coupled seepage–stress analyses. In addition, the rainfall characteristics (intensity, duration, antecedent weather condition, etc.) of a landslide slope also have significant effects on the slope response. Ng and Shi (1998), Kristo et al. (2017) and Rosly et al. (2023) pointed out that the pore pressure increases quickly, and the FOS decreases drastically when infiltration rate comes close to soil permeability. Rahimi et al. (2011) performed parametric seepage–stability analyses to show that different antecedent rainfall patterns markedly alter pore-water pressure evolution and can significantly reduce the factor of safety, especially in low-conductivity slopes. Moreover, Wu et al. (2015b) developed a physical model test that indicated well-developed fissures could accelerate concentrated infiltration, which produced a local and fully seepage zone and lead to a preassembled slide plane. With the development of computational geomechanics and information, numerical modeling has played an important part in elucidating rainfall-seepage-failure processes. For example, Cai and Ugai (2004), used the strength reduction finite element method to simulate the process of slope stability evolution induced by rainfall infiltration. Zhang et al., 2005 studied the evolution of the factor of safety under rainfall in a time step with finite element and strength reduction method. Although these methods could accurately simulate the pore pressure and stress changes during rainfall, they are incapable of reproducing significant deformations after the failure of the slope because of mesh distortion. This shortcoming can be remedied by introducing particle-based methods to rainfall triggered landslides research. Wang et al. (2024), based on a GPU accelerated material point method, proposed a hydro-mechanical coupled computation scheme that can realize the surface runoff model and the subsurface seepage model which can continuously simulate slope deformation and failure at every instant of rainfall.

Beyond general rainfall behavior, some research has shown that engineering structures and anthropogenic surfaces can significantly modify rainfall infiltration pathways. Surface covers, impervious facilities, and engineered drainage systems (Basher and Ross, 2001; Jordán and Martínez-Zavala, 2008; Kalantari et al., 2014) can alter hydrological boundary conditions and lead to localized seepage. These findings highlight the importance of accurately modeling spatially variable rainfall inputs when assessing slope stability under realistic conditions. In mountainous photovoltaic (PV) power stations, this type of human-induced modification to rainfall distribution becomes particularly evident. In recent years, large-scale PV power stations have been widely constructed in mountainous regions. While providing clean energy, they introduce new geohazard challenges because PV panels create spatially non-uniform rainfall distribution through shading and runoff concentration. Such heterogeneity significantly alters infiltration patterns and challenges the applicability of conventional rainfall-induced slope stability methods. However, limited studies have specifically investigated the hydro-mechanical impacts of PV-induced rainfall redistribution.

However, most existing studies consider a spatially uniform rainfall distribution and are therefore unable to properly account for localized non-uniform infiltration effects due to rainfall redistribution under PV facilities. Systematic efforts in this direction are lacking. Overall, although many advances have been made regarding the unsaturated seepage behavior and the effect of rainfall, there still exist three important issues that deserve further studies: (Ⅰ) Most existing hydro-mechanical analyses assume spatially uniform rainfall, which fails to capture the true hydrological conditions experienced on PV-covered slopes. (Ⅱ) Under complex topography and non-uniform boundary conditions, the temporal characteristics of coupling mechanism between infiltration, pore pressure and shear strength still lack a precise description. (Ⅲ) The change in soil parameters and partial erosion during the unsaturated-saturated change was rarely taken into consideration in overall stability analysis.

Based on spatially non-uniformity rainfall on the slope and influence of rainfall frequency in sloping photovoltaic regions in mountainous areas, this paper proposed a hydro-mechanical coupling model. Using FLAC3D numerical simulation and theory analysis, this paper systematically explores the evolution of pore water pressure generated by rainfall infiltration and its effect on slope stability. For the first time, our work calculates the non-uniform rainfall effects caused by PV panel layouts on the slope stability analysis, including the corresponding hydro seepage responses and failure processes. The results provide an analytic insight and practical reference to slope protection design and slope-warning systems of mountainous PV engineering projects.

2 Theoretical framework and research methodology

2.1 Unsaturated-saturated seepage modeling

Rainfall infiltration is one of the essential mechanisms in the water cycle (Fleckenstein et al., 2010; Sophocleous, 2002). The water infiltration from rainfall to slope surface gradually rises to the volumetric water content of soil. Once the near-surface layer of soil reaches a saturated state, part of infiltrated water will become surface runoff. Changes in the slope geometry or surface cover such as shading of the PV panels introduce pronounced spatial heterogeneity of infiltration, resulting in a non-homogeneous seepage field (Figure 1).

Figure 1

Figure 1. Schematic diagram of rainfall-infiltration-runoff process on slope.

Because of the dominant influence of liquid-phase flow on slope stability, we adopt a single-phase seepage model to simulate the unsaturated–saturated flow process of soil and neglect air pressure and air–water interactions for simplicity, since their influence on short-term rainfall-induced stability is smaller than that of liquid-phase flow. This single-phase assumption has been widely used and validated in analyses of rainfall-induced landslides (Lu and Likos, 2004; Ng and Shi, 1998; Rahimi et al., 2011).

2.1.1 Governing equations and constitutive relationships

From Darcy’s law and the continuity equation, the governing equation of unsaturated seepage can be obtained, giving rise to the famous Richards equation (Richards, 1931). This equation is the fundamental equation of the transient unsaturated-saturated seepage flow in porous media. The governing equation for unsaturated seepage is given in Equation 1.

$\frac{\partial \theta }{\partial t}=\nabla ·\left[K\left(\theta \right)\left(\nabla \psi +i_z\right)\right](1)$

where θ denotes the volumetric water content, t represents time, K is the unsaturated hydraulic conductivity, ψ denotes the pressure head, and i_z is the unit vector in the gravitational direction.

To more accurately capture the coupling between the volumetric water content and pressure head during rainfall infiltration, the mixed form of the Richards equation is used to simultaneous moisture continuity and hydraulic potential variations such as:

$\frac{\partial \theta }{\partial t}=\nabla ·\left[K\left(\theta \right)\left(\nabla \psi +i_z\right)\right]=C\left(\psi \right)\frac{\partial \psi }{\partial t}(2)$

When K(θ) is the unsaturated hydraulic conductivity, and C(ψ) represents the specific moisture capacity, or the unit increases in the volumetric water content corresponding to a variation of 1 unit of pressure head. It is a key index in Equation 2 that expresses the transition from unsaturated to saturated infiltration due to rainfall infiltration. The functional form of this equation is often given by Equation 3:

$C\left(\psi \right)=\left({\theta }_s-{\theta }_r\right)\alpha n{\left[1+\left(\alpha \left|\psi \right|\right)^n\right]}^{-\left(1+1/n\right)}(3)$

2.1.2 Soil-water characteristic and hydraulic conductivity models

Hydraulic conductivity in unsaturated soils has both a volumetric water content and matric suction effects. It is therefore characterized by using the VG model (Mualem, 1976; Van Genuchten, 1980) of hydraulic functions, as shown in Equation 4. It is an empirical functional form: modelling the soil–water characteristic curve (SWCC) as follows, enabling a smooth relationship between pressure head and degree of saturation:

$S_e=\frac{\theta \left(\psi \right)-{\theta }_r}{{\theta }_s-{\theta }_r}={\left[\frac{1}{1+{\left(\alpha \left|\psi \right|\right)}^n}\right]}^{1-\frac{1}{n}}(4)$

Modelling of relation of unsaturated hydraulic conductivity to relative saturation was usually expressed in the form of Mualem-modified VG model, a simple hydraulic law that reflects the pore structure effect and effect of conductivity change with saturation. The hydraulic conductivity function corresponding to the modified VG model is given in Equation 5. The modified Mualem model is usually used to give a better precision of infiltration prediction in unsaturated soils:

$K\left(\theta \right)=K_s\times {S_e}^{1/2}{\left[1-{\left(1-{S_e}^{1/n}\right)}^{n/\left(n-1\right)}\right]}^2(5)$

where S_e denotes the relative saturation, θ_s is the saturated water content, θ_r is the residual water content, Ks is the saturated hydraulic conductivity, and n, α are empirical coefficients. Compared with the traditional VG model, the modified expression has better stability in the area with high suction and is widely used for seepage analysis in an unsaturated-saturated transition.

Governing equation for unsaturated seepage and the soil-water characteristic model lay the theoretical basis of simulating unsaturated conditions and pore-pressure field spatially non-uniform generated by the effect of PV shading. This theory provides more precise representation of the hydrological response for the uneven rainfall over the whole slope by coupling of hydrological and geo-mechanical processes.

Rainfall-induced slope instability is a transient and dynamic process. As soon as rainfall starts, infiltrating rainwater gradually fills soil pores and infiltrates downward with increasing volumetric water content and bulk density (Zhang M. et al., 2025). As the matric suction in the unsaturated zone decreases, the effective stress generated by the soil’s self-weight also diminishes, leading to a reduction in shear strength and a subsequent decrease in slope stability (Song et al., 2021). As more rainwater infiltrates, the dry soil layers become saturated, and pore-water pressure builds up gradually, and eventually develops to positive pore pressure at the surface. The total stress (Zienkiewicz et al., 1990) equilibrium equation is expressed as follows.

$\frac{\partial {\sigma }_{ij}}{\partial x_j}+{\rho }_s{g}_i=0(6)$

The bulk density of unsaturated soil becomes denser as rainfall seeps into soil. The equation for calculating the bulk density is shown as follows:

${\rho }_s={\rho }_d+ns{\rho }_w(7)$

where ρ_s is the saturated density, n is porosity, and s is the saturation degree (the saturation degree can be calculated through solid density ρ_d and fluid density ρ_w, and the effective stress is defined as follows:

${{\sigma }_{ij}}^{\prime }={\sigma }_{ij}+\alpha p{\delta }_{ij}(8)$

Substituting Equations 7, 8 into Equation 6, then the equilibrium equation for soil structure under saturated conditions is shown as in Equation 9.

$\frac{\partial {{\sigma }_{ij}}^{\prime }}{\partial x_j}+{\rho }_s{g}_i-\alpha \left(1-n\right)\frac{\partial p}{\partial x_i}-\eta {\gamma }_w\frac{\partial \varphi }{\partial x_i}=0(9)$

where α is the Biot coefficient, n denotes the porosity, p represents the pore-water pressure, γ_w is the unit weight of water, Φ denotes the hydraulic head, and g_i is the gravitational acceleration.

2.1.3 Shear strength under unsaturated conditions

Because matric suction exerts a significant influence on soil strength, rainfall-induced changes in pore-water pressure directly cause the redistribution of effective stress and shear strength. To quantitatively describe this relationship, Fredlund DG. and Rahardjo H. (1993) proposed the dual stress state variable model based on the principle of effective stress, which is expressed in Equation 10 as follows:

${\sigma }^{\prime }=\left(\sigma -u_a\right)+\chi \left(u_a-u_w\right)(10)$

The corresponding shear strength under unsaturated conditions is given by Equation 11.

${\tau }_f=c^{\prime }+\left(\sigma -u_a\right)\mathit{tan}{\phi }^{\prime }+\left(u_a-u_w\right)\mathit{tan}{\phi }_b(11)$

In this equation, σ′ represents the effective stress, σ denotes the total stress, u_a and u_w are the gas-phase and water-phase pore pressures, respectively, χ is the suction stress distribution coefficient (which is related to the volumetric water content), c' is the effective cohesion, φ′ is the effective internal friction angle, φ_b is the friction angle associated with matric suction, and τ_f is the shear strength at failure. For simplification, this study neglects air-phase effects, assuming matric suction to be equal to pore-water pressure. When the soil is in an unsaturated state, the strength generated by matric suction can be included in the cohesion. The total cohesion of the soil can be expressed in Equation 12 as:

$c^{″}=c^{\prime }+\left(u_a-u_w\right)\mathit{tan}{\phi }_b(12)$

Based on the SWCC, the matric suction corresponding to different volumetric water contents can be directly determined. Consequently, the cohesion can be adjusted accordingly to reflect the changes in shear strength of the soil under unsaturated conditions. When the soil becomes saturated, the total normal stress equals the sum of the effective stress borne by the soil skeleton and the pore pressure carried by pore water, following Terzaghi’s effective stress principle (Terzaghi et al., 1996), the total normal stress under saturated conditions can be expressed by Equation 13:

$\sigma ={\sigma }^{\prime }+u_w(13)$

The pore-water pressure uw can be computed using Equation 14:

$u_w={\rho }_{w}gh(14)$

where σ is the total normal stress, σ′ is the effective normal stress, u_w is the pore-water pressure, and h is the positive hydraulic head. When the soil nodes are saturated, the effective stress is calculated using the above formula.

2.2 Methods for stability analysis of slopes under rainfall conditions

Strength reduction method (SRM) is the most widely used numerical method to calculate FOS of slopes. Its basic idea is to uniformly reduce the shear strength index (cohesion and friction angle) of the soil until it reaches equilibrium or breakdown state on the limit surface of the slope. The reduction factor at this point is defined as the safety factor. Based on the reduced effective internal friction angle (φ_F) and effective cohesion (c_F), the reduced shear strength can be expressed in Equations 15–17 as:

$c_f=\frac{c}{R_f}(15)$

$\mathit{tan}{\phi }_f=\frac{\mathit{tan}\phi }{R_f}(16)$

${\tau }_f=c_f+\sigma \mathit{tan}{\phi }_f=\frac{c}{R_f}\sigma \mathit{tan}\left[{\mathit{tan}}^{-1}\left(\frac{\mathit{tan}\phi }{R_f}\right)\right](17)$

When performing strength reduction under unsaturated seepage conditions, the effective stress principle must be integrated to account for the real-time effects of pore-water pressure. During the reduction process, the shear strength of each element is calculated based on its current effective stress state. Therefore, the increase in pore pressure caused by rainfall will directly result in a decrease in shear strength and a corresponding reduction in slope stability after strength reduction.

2.3 Custom numerical methods for coupling seepage and mechanics in slope analysis

This study employs the consolidation–seepage module in FLAC3D and integrates custom-developed FISH scripts to achieve coupled analysis of unsaturated-saturated seepage and mechanical responses. The built-in flow–solid coupling module in the software primarily targets saturated soils and is unable to directly simulate the dynamic evolution of matric suction and unsaturated zones induced by non-uniform rainfall infiltration. To address this, the study extends the module’s functionality by developing custom FISH scripts, enabling the calculation of unsaturated soil slope seepage. Figure 2 illustrates a schematic diagram of this coupled algorithm.

Figure 2

Figure 2. Coupling framework of rainfall-hydrology-geomechanics model.

2.3.1 Seepage calculation module for unsaturated-saturated flow

In the software seepage analysis, the fluid flow follows Darcy’s law, as shown in Equation 18.

$q_i=-k_1\widehat{k}\left(s\right){\left[p-{\rho }_f\gamma g_j\right]}_j(18)$

where q_i represents the flow vector in the specified direction, p is the pore-water pressure, and k_ij is the conductivity tensor. And $\widehat{k}$ (s) denotes the relative permeability, defined as s²(3-2s), where s is the saturation degree, ρ_f is the fluid density, and g_j (j = 1, 2, 3) represents the three components of gravitational acceleration. During seepage calculations, FLAC3D executes either saturated or unsaturated flow modules independently. In saturated seepage, the software sets the saturation degree at each node to 1, and the flow coefficient is forced to the saturated permeability coefficient. The formula for updating the pore-water pressure at the next step is, as shown in Equation 19:

$\frac{1}{M}\frac{\partial p}{\partial t}=\left(-q_{i,i}+q_v\right)-\alpha \frac{\partial \epsilon }{\partial t}(19)$

where M is the Biot modulus (Pa), and n represents the porosity. When the compressibility of the soil is not considered, α is set to 1. In the case of unsaturated seepage, the software sets the pore-water pressure at each node to zero, and the update of saturation is computed using Equation 20:

$\frac{s}{n}\frac{\partial s}{\partial t}=\frac{1}{s}\left(-q_{i,i}+q_v\right)-\alpha \frac{\partial \epsilon }{\partial t}(20)$

where s represents the degree of saturation, n is the porosity, and ε denotes the volumetric strain.

2.3.2 Calculation of mechanical parameter modifications

When the soil becomes saturated, the total normal stress consists of the effective stress carried by the soil skeleton and the pore-water pressure borne by the pore fluid. The strength and deformation of the soil are governed by the effective stress, which determines its mechanical behavior. This concept forms the basis of the effective stress principle for saturated soils proposed by Terzaghi et al. (Terzaghi et al., 1996). Once the soil becomes saturated, the effective stress can be evaluated using Equation 21:

$\sigma ={\sigma }^{\prime }+u_w,u_w={\rho }_{w}gh(21)$

where σ is the total normal stress, σ′ is the effective normal stress, u_w is the pore-water pressure and h is the hydraulic head. Based on the Richards equation and SWCC characteristics, a FISH function was generated to modify the saturation degree and unsaturated permeability coefficient of each element at each time step. Specific steps for this procedure are described as follows. At each sub step in solving the flow–solid coupling equations, for each element read its current pore-water pressure and get its saturation degree and relative permeability from the saturation–conduction degree relative permeability curve given by the Van Genuchten model. Following this, the permeability coefficient of this element will be dynamically adjusted. At the same time, the calculated degree of saturation is assigned to nodes to define the surface saturation state. This module can make model changes in soil permeability properties due to rainfall infiltration, to achieve unsaturated seepage calculations using FLAC3D.

2.3.3 Rainfall boundary condition application and transition

To simulate spatially non-uniform rainfall infiltration across the slope surface, FISH scripting is employed to assign distinct infiltration boundary conditions to different slope zones. In the initial state, the entire slope surface was defined as an infiltration zone. Then, according to the horizontal projection range of the PV panels, the corresponding slope sections were reclassified as shaded zones, as illustrated in Figure 3.

Figure 3

Figure 3. Water flow paths and rain-shielded areas under photovoltaic panels on the slope.

First, the model identifies the element surfaces belonging to infiltration and shaded zones. For infiltration-zone elements, a rainfall infiltration flux is applied, while for shaded-zone elements, a lower infiltration flux is assigned to represent the limited rainfall entering through panel gaps or lateral splashing. In addition, a dynamic boundary adjustment mechanism is incorporated into the unified rainfall–seepage iteration process for both zones. If a positive pore-water pressure develops on any infiltration boundary element, the infiltration process at that location is temporarily halted, and the boundary condition is switched to a zero pore-pressure boundary. Such treatment causes surface runoff to take place or for there to be ponding in soils, with excess water being directed to wetting around neighboring unsaturated regions for further infiltration, and once the local pore pressure reaches a negative value (the ponded water has dissipated), the infiltration flux boundary is switched on again.

By means of this dynamic boundary control strategy, any slope elements can be assigned either a specified-flux or a specified-head boundary condition during the rainfall process. A series of logical scripts enables the spatial temporal updating of rainfall infiltration boundary conditions across the entire slope surface, providing a more realistic representation of the evolution of seepage boundaries under non-uniform rainfall conditions.

2.3.4 Strength reduction and stability criterion module

Fredlund and Rahardjo (1993b) developed a shear strength equation for unsaturated soils using a dual stress state variable approach. The shear strength equation is expressed in Equation 22 as follows:

$\tau =c^{\prime }+{\left(\sigma -u_a\right)}_f\mathit{tan}{\phi }^{\prime }+{\left(u_a-u_w\right)}_f\mathit{tan}{\phi }^b(22)$

where c′ and φ′ represent the effective cohesion and effective internal friction angle of the soil, respectively. σ denotes the total normal stress acting on the failure plane at the point of failure. ua and u_w correspond to the pore-air pressure (assumed to be atmospheric pressure in this study) and the pore-water pressure on the failure plane, respectively. (σ−u_a)_f and (u_a−u_w)_f denote the net normal stress and the matric suction at failure, respectively, while φ^b represents the angle indicating the rate of change of shear strength with matric suction. The shear strength in the unsaturated zone is corrected by introducing matric suction and its associated friction angle, which yields a generalized cohesion term, expressed in Equation 23.

$c^{*}=c^{\prime }+{\left(u_a-u_w\right)}_f\mathit{tan}{\phi }^b(23)$

To analyze the stability of the slope at different rainfall stages, the “averaging” method was used. The rainfall process was divided into several periods. At the end of each period, the stress–pore pressure field of the model was extracted and taken as the initial condition of the following strength reduction analysis. Using a FISH script of batch processing mode, the program automatically saved the model state after one stage of rainfall–seepage calculation is completed, and removes the flow–solid coupling function, then switches to strength reduction mode to calculate its FOS. After the stability analysis is completed, the model reactivates the flow–solid coupling stage, loads the initial pore pressure value in the next stage again, and continues the rainfall infiltration simulation stage. This cycling process goes on until all the period of rainfall is simulated, and this completes the whole period of seepage–strength reduction analyses and makes the hydraulic and mechanical response of the slope coupled with each other throughout the whole period of rainfall process. Summarizing the above discussion, our secondary development produces a coupled numerical calculation model that considers the influence of non-uniform rainfall on slope instability, whose method solves the variation field of seepage flow rate and the stress–deformation field of soil at each rain time step simultaneously. By combining explicit boundary condition treatments with strength-reduction analysis, the method also enables full tracking of the evolution of the internal seepage–stress field and gradual evolution of the slope stability in rainfall events. Our approach differs from previous studies in that it utilizes the built-in flow–solid coupling capability of the software, together with newly developed modules, to model complex rain scenarios within a finite-difference framework. Consequently, the model can also include the rainfall redistribution effect of photovoltaic panels without introducing another individualized computational scheme.

3 Overview of the research area

Based on the algorithm put forward in this paper, we research the stability and failure of a slope inside a PV station in a west mountain area under uniform rainfall and spatially non-uniform rainfall conditions respectively. This slope shows a common trend from southwest to northeast and with elevation gradually increasing from southwest to northeast between 1386 m and 1434 m in position. Conversely, the overall terrain is relatively flat, while several local spots have steep slopes and obvious cut depressions of gully.

As shown in Figure 4, the typical geological section of slope consists of three strata. The surface layer is Quaternary eluvial deposits (Q₄^el) which are mostly silty clay intercalated with gravel. This layer is loose in texture, relatively large pores are developed, therefore its permeability is medium. Below is a clay transition layer functioning as a dense and impermeable transition zone between the loose overlying layers and the underlying bedrock. This transitional layer weathers to a non-resistance zone once water seeps into it. The lower bed is a strongly weathered bed of basalt (P₂β), with intense weathering and low strength.

Figure 4

Figure 4. Geological cross-section along line F–F′ showing stratigraphy of the PV-site slope.

On the rainfall regime of the PV power station area, strong seasonality is a distinguishing feature, featuring torrential and intense rain events lasting for only a few hours most of the time in summer. The hydrogeological properties in this PV power station site are rather plain, mainly because groundwater recharge is replenished from rainfall in the atmosphere. None of the boreholes took up groundwater, thus suggesting that the water table lies at a great depth. The annual rainfall totals from our study area are given in Figure 5, which we can use to define annual rainfall boundary conditions in the following seepage–stability analysis.

Figure 5

Figure 5. Annual rainfall and duration curves of the PV-site slope used for defining rainfall boundary conditions.

Field investigation showed that some basic drainage ditches were located under the back of the PV panels to reduce the runoff erosion at the slope toe by thin sheets of water. However, the whole slope could still be damaged by local slides or over-turn of very steep rainfall. The arrangement of the PV panels also leads to rather strong rainfall redistribution on the slope ground, forming alternate wet and dry zones. Soil beneath the panels remains relatively dry because of shading. But in the exposed inter-panel soil, soil moisture becomes increased markedly during rainfall. After rainfall events, visible deformations were formed in the exposed inter-panel areas, large shallow scale collapses happened across the slope (Figure 6).

Figure 6

Figure 6. Schematic of the investigated PV slope showing drainage layout and rainfall-induced surface deformation zones.

These observations show clearly that photovoltaic panels affect rain infiltration and slope movement significantly. So, a numerical model is needed to understand the slope stability under rainfall–shading effect induced by PV panel arrays. From the geological and geomorphological characteristics of the PV site and its field results, the 3D numerical model of the slope was developed (Figure 7a). The modelled slope is 45.98 m and height 36.46 m. To ensure the accuracy as well as computational efficiency of the modelling, the thickness of the model perpendicularly to the slope is set to 2 m. The mesh has a total of 69,790 nodes and 63,665 elements. Boundary conditions were applied with zero displacement conditions across the front, back and bottom boundary in the vertical direction, and the lateral boundary to horizontal direction; zero displacement constraints were applied across both lateral boundaries in horizontal direction. The soil was initially assumed to be in an unsaturated state, with the pore-water pressure field defined to represent in situ hydrostatic conditions prior to rainfall infiltration. The Mohr–Coulomb constitutive model was employed to simulate the mechanical behavior of slope materials. All model parameters were determined through laboratory tests and field investigations, as summarized in Table 1. Following parameter ranges reported in previous studies (Bryant, 2003; Garcia-Bengochea et al., 1979), the saturated hydraulic conductivities of the silty clay and clay layers were set to 5 × 10⁻⁷ m/s and 5 × 10⁻⁸ m/s, respectively, whereas that of the strongly weathered basalt was assigned as 1 × 10⁻⁶ m/s. Although hydraulic properties of geomaterials can vary spatially and may evolve due to wetting–drying cycles (Nguyen et al., 2022; Jing et al., 2022; Liu et al., 2024), the three formations considered in this study all exhibit fine-grained structure and similar pore size distribution, which is consistent with the borehole description and laboratory test results. Therefore, in the present analysis, we assign a specific layer of saturated hydraulic conductivity but assume a common set of VG parameters for the SWCC and the relative permeability function. In addition, the duration of the rainfall events analyzed in this paper is short, so hydraulic hysteresis and long-term microstructural evolution are expected to be limited. Under these conditions, the use of a unified VG model for all formations is a reasonable simplification (Lu et al., 2012a).

Figure 7

Figure 7. Three-dimensional numerical model of the PV-site slope. (a) The scheme of the geometric model (b) Comparison of uniform and non-uniform rainfall infiltration boundary conditions.

Table 1

Table 1. Material parameters of the slope.

Because the primary hydrological effects of the PV support system arise from panel shading and runoff concentration, non-uniform rainfall distribution on the slope surface was simulated by adjusting the infiltration boundary conditions of corresponding surface elements. According to field observations, alternating shaded and unshaded zones were delineated along the model slope surface. Each PV panel row covers a projected length of approximately 4.4 m, while the exposed interval between adjacent rows is about 1 m (Figure 7b). The selected rainfall intensity of 1.157 × 10⁻⁶ m/s (uniform rainfall) and 3.47 × 10⁻⁶ m/s (non-uniform rainfall) represents the extreme 24-h rainfall recorded in the PV station, corresponding to a return period of 20–50 years. These values are consistent with hydrological records from 2018 to 2025 and are used in regional slope design. A 24-h continuous rainfall event with an average intensity of approximately 100 mm/day was applied to represent short-duration heavy rainfall typical of summer storms in the study area. To implement this spatial variability, non-uniform rainfall boundary conditions were implemented by assigning different infiltration fluxes to the exposed and PV-shaded areas. The infiltration zone was prescribed to a flux of 3.47 × 10⁻⁶ m/s, corresponding to 80% of the incident rainfall, while the shaded area received only 0.10 × 10⁻⁶ m/s, representing the residual dripping and micro-runoff that bypasses panel edges. Following the procedure described in Section 2.3, the corresponding infiltration flux was applied to infiltration zones, while shaded zones were assigned a substantially smaller flux, reflecting limited indirect precipitation. Such configuration results in the overall balance of rainfall over the slope and replicates concentrations of runoff and spatial variation in infiltration that results from specific PV panel arrangement in space.

The non-uniform rainfall boundary is applied by a FISH routine that assigns time-dependent flux values to surface nodes based on the shading–infiltration mapping derived from the PV panel geometry. The routine only updates the boundary flux at each timestep, while the coupled flow–mechanical equations are fully solved by FLAC3D’s built-in algorithm. Under uniform rainfall input, the model reproduces pore-pressure evolution trends consistent with results reported by (Wang et al., 2024) Under no-rainfall conditions, the system maintains hydrostatic equilibrium, confirming that the FISH routine introduces no artificial flux. These steps confirm that the numerical framework is stable and reliable for evaluating the effects of PV-induced non-uniform infiltration.

4 Analysis of slope stability under uniform rainfall condition

4.1 Characteristics of rainfall infiltration

Figure 8 shows the infiltration depth evolution for the slope under uniform rainfall condition. The infiltration of the slope goes through three step-up processes: rapid stage, transition stage and quasi-steady stage, representing different features of slope hydraulic change. At the initial period (0–5 h), the infiltration depth increased very fast between 0 and 1.67 m within the initial 0–5 h, which has an average wetting-front velocity of 0.334 m/h. The rapid wetting speed was attributable to the accumulative action of initial matric suction field and the high permeability of silty clay in near surface. When the surface soil almost reached saturated state, its hydraulic conductivity increased with the degree of saturation, and a rapid infiltration phase appeared at the beginning of rainfall. Within 5 h–15 h, the depth of infiltration increased only slightly (by 0.54 m), the transition into slower infiltration regime. The sudden fall in the permeability coefficient and the difference in stiffness between the two layers caused perched water accumulation—infiltration here was jointly limited by low permeability and restricted drain ability. At 15 h, the infiltration was completed to a quasi-steady state (from 15 to 24 h). Rainwater infiltration proceeded slowly, mostly in the gravity-flow and diffusion-controlled regime. A high saturation zone developed near the slope surface and moved downslope and reached towards the toe, where it finally intersected the weak basal layer.

Figure 8

Figure 8. Evolution of infiltration depth in the slope under uniform rainfall conditions.

4.2 Change in stability

Figure 9 presents the temporal evolution of the FOS of the slope under uniform rainfall. As rainfall duration increased, the FOS exhibited three distinct stages of variation: a rapid-decline stage (0–5 h), a steady-decline stage (5–15 h), and a pre-failure stage (15–24 h). In the rapid-decline stage, the FOS dropped sharply from 2.45 to 2.27, with an average reduction rate of approximately 0.036 h^-1. Because the rainfall infiltration rate exceeded the transmissivity of the upper soil layers, excess water accumulated near the slope surface, forming a high-moisture zone. The matric suction in both the silty clay and clay layers dissipated rapidly, resulting in a progressive reduction in equivalent cohesion and internal friction angle. During the steady-decline stage (5–15 h), the FOS further decreased to 2.13 and the decline rate slowed. The advancing wetting front advanced to the upper edge of the clay transition region and gave rise to a perched water zone with accumulated local excess pore-water pressure forming in it. The shallow slope section was approaching saturation, and the happening rainfall travelled only at slow speed for a downward migration of the wetting front. Finally, at the pre-failure stage 15–24 h, the rate of FOS reduction increased again. As the analysis progressed and the fluids infiltrated, an accumulation of positive pore-water pressures developed in the neighborhoods of the slope toe, which coalesced with the downdip (sliding face) oriented shear band developing in the upper layers and caused a rapid degradation in slope stability.

Figure 9

Figure 9. Temporal evolution of the FOS of the slope under uniform rainfall conditions.

Figure 10 shows the plastic strain contours of the slope for different times. The plastic shear band first occurred in the upper silty clay layer, forming a thin strip along the slope surface. With the continuous rainfall infiltration, the shear band propagated down the slope, and the magnitude of plastic shear bands grew as they became wider. By 15 h, the plastic reached the interface between the silty clay layer and the clay layer, which formed a relatively continuous weak zone; at 24 h, the plastic zone connected with the weak zone near the foot of the slope, forming a basically continuous potential sliding surface, reaching from the crest to foot of slope; at this time the FOS decreased to 1.95, and the slope is critically in the critical state of stability. Just as this spatial evolution of the shear zone matches the trend in the temporal FOS, this confirms progressive suction loss and local reductions in shear strength as the main contributors to rainfall-induced slope failures.

Figure 10

Figure 10. Progressive development of maximum shear-strain rate and potential slip surface under uniform rainfall.

Figure 11 illustrates the displacement evolution of the slope under various rainfall durations, revealing a typical shallow translational deformation pattern. Deformation was mainly concentrated in the upper shallow zone near the slope crest and in the toe region. During the initial stage (1–5 h), corresponding to rapid suction dissipation, rainfall infiltration was confined to the upper 1 m of the silty clay layer. Matric suction near the surface decreased sharply, producing small downslope displacements and slight outward extrusion at the slope toe. During the second stage (5–15 h), the wetting front migrated downward, and the displacement peak gradually shifted from the mid-slope toward the toe. By 24 h, the wetting front penetrated the base of the clay layer, reaching the interface with the strongly weathered basalt. Integrating the characteristics of all three stages, the reduction in matric suction and the buildup of pore-water pressure emerge as the dominant mechanisms controlling displacement evolution. Moreover, the increasing volumetric water content increased the unit weight of the soil, further boosting downslope sliding tendencies and increased the slope susceptibility to failure during, or soon after, heavy rain. Strength–loading path interactions were largest within the two upper layers, where loss of suction and the pore-pressure increase reduced shear resistance together. However, the underlying strongly weathered basalt was a competent bearing layer that provided support for a foundation but did not participate directly in failure. Overall, the shallow soil layer was highly sensitive to rainfall infiltration, whereas the bedrock was more affected by structural discontinuities than by hydraulic effects.

Figure 11

Figure 11. Evolution of cumulative displacement fields under uniform rainfall infiltration.

5 Comparative analysis of slope stability under non-uniform rainfall conditions

5.1 Characteristics of rainfall infiltration

Figure 12 illustrates the change of saturation field with non-uniform rainfall, in which we find distinct stripes and up-and-down fluctuation. In the infiltration domains, the matric suction in the surface silty clay diminishes first, and hydraulic conductivity will enlarge suddenly with increasing saturation. Consequently, the wetting front advances downward in a preferential, banded manner, leading to early developments of saturated zones everywhere. The contrast is that on its neighboring, shaded patches with a weaker intensity of rainfall, soil suction is still relatively high, causing saturation increase to be delayed. This contrast forms a distinct lateral hydrostatic potential difference between the infiltration zone and shaded ones. As a result, the saturation distribution is wavy and alternating, peaking at the zones of infiltration where local runoff cones quickly accelerate early saturation. As infiltration flows on, neighboring areas of high saturation gradually connect up at the mid-slope area. By 8 h, the slope surface becomes nearly saturated; however, the infiltration front is spatially still rather uneven, showing oscillating progress in the direction of the slope. The wetting front below the infiltration zones is clearly deeper, while it is lagging behind below shaded zones, this revealing a clear spatial synchronicity of the infiltration.

Figure 12

Figure 12. Temporal evolution of the saturation field in the slope under non-uniform rainfall conditions.

Under uniform rainfall, the saturation boundary is nearly planar, and the wetting front propagates parallel to the slope surface. But after 9 h, differences in infiltration due to non-uniform rainfall vanish gradually, because the longer rainfall and redistribution of moisture eventually partially balance out the difference caused earlier. Overall, non-uniform rainfall changes the saturating evolution pattern from a nearly translational advancing pattern to a banded and merging pattern laterally. Localized saturation starts earlier and progressively converged at the slope toe. Earlier stages of rainfall: soil below the shaded area also eventually become saturated in the later stages of rainfall as lateral diffusive moisture flows into saturated adjacent infiltration zones. Those pore-water pressure observations support this change: sensors in a main infiltration region measure rapid increase and prompt stabilization; sensors in the shaded region respond later but slowly catch up or indeed go above them if they occur in rainy conditions. Non-uniform rainfall leads to a pronounced hydrological heterogeneity. Spatial and temporal inconsistencies of saturation and pore-pressure developments mean spatially different decreases in shear strength, with possible instability manifestations at early times in localization.

Figure 13 shows the temporal evolution of the infiltration depth distribution along the slope for non-uniform rainfall conditions. It can be seen that the infiltration process is much more scattered with temporal evolution, especially at the initial stage of the rainfall condition. Within the first 1–5 h, the wetting front under the infiltration zones developed quickly downwards, with a significantly higher infiltration rate than in the shaded zones. This process resulted in a highlighted wave-shaped infiltration profile along the slope. As can be seen from the enlarged image above, at 1 h the infiltrations depth reached 1.0 m approximately in the infiltration area and reached beyond 1.6 m at 5 h, but depths in the nearby shaded area are less than 1.0 m. As the rainfall continued, the wet fronts in the slope continued to grow and link to each other, the spatial difference between them was further reduced. By 24 h, the overall infiltration front generally settled down, and the amplitude of the shape-wave front decreased much. By contrast, under uniform rain intensity, the spatial distribution of infiltration depth was smooth and monotonic, and the wetting front progressed steady and parallel to the underlying slope surface. Non-uniform rain was able to transform the behaviour from a smooth and monotonic distribution of infiltration depth into a fluctuating and stratified manner and with a peculiar lateral asynchronous fluctuation in subsurface waters. These hydrological heterogeneities formed for two reasons. First, in the infiltration zones where the rainfall was higher, the matric suction dissipated fast, increased hydraulic conductivity and the wetting front propagated fast. Second, in the shaded areas where the rainfall inputs were limited, the drying was slow so the suction dissipated delayed, and infiltration was slow. The contrast also created a pronounced moisture potential gradient between the two zones, leading to lateral seepage recharge from wetter to drier zones. With continued rain, the influence of gravity flow and infiltration from the side of components reduce this difference gradually. However, at the toe zone and downstream at major flow outlets, the front of the infiltration front is deeper and pore-water pressure is much higher, and the local hydraulic concentration has yet to dissipate.

Figure 13

Figure 13. Distribution of infiltration depth at different rainfall moments.

The non-uniform rainfall pattern disrupted the continuity of infiltration depth along the slope and promoted the premature formation of high-saturation and high pore-pressure zones within runoff convergence regions. This pronounced spatial variability in infiltration serves as a key hydrological control factor governing the subsequent degradation of slope stability and the early development of continuous shear bands. These hydrological contrasts establish the direct physical foundation for the following section, which analyses the evolution of the factor of safety and displacement response under non-uniform rainfall conditions.

5.2 Stability evolution and instability mechanism analysis

To study the surface deformation under non-uniform rainfall, we set three monitoring points at the toe (P1, x = 8.08 m), mid-slope (P2, x = 28.84 m) and crest (P3, x = 34.34 m). The corresponding displacement–time curves are shown in Figure 14.

Figure 14

Figure 14. Evolution curve of displacement with time at different monitoring points.

Results indicated that displacements at the three points all share the same three-stage development pattern. When the slope enters the onset and localized wetting stages, the slope mainly represents the surface wetting and a small amount of surface deformation from local infiltration. At this subsequent period of rapid growth, when the positive pore-water pressure in the infiltration zones keeps growing and spreading to below the gradient section, the displacements at the three monitoring points increase abruptly and almost simultaneously. Meanwhile, a lateral hydraulic gradient develops beneath the infiltration zones, accompanied by differential deformation, causing the seepage lines to converge progressively toward the slope toe. Meanwhile, a lateral hydraulic gradient occurs beneath the infiltration zones, and differential deformations also occur, so that seepage lines converge to slope toe progressively. After about 17 h, displacement rate decreases sharply, the deformation rate drops from rapid plastic motion to a quasi-steady state. Eventually, the cumulative surface displacements at P1, P2 and P3 are found to reach about 6.64 mm, 5.46 mm and 5.41 mm, respectively. All these reveal that non-uniform rain does not increase the total amount of infiltration, but rather further increases the spatial heterogeneity of both the seepage and the deformation responses. As a result of this uneven hydromechanical behavior, there is locally coarse early slope softening, in particular in the area of the slope toe, and a precursory pattern of instability in terms of concentrated deformation and hydraulic response in the lower part of the slope.

Figure 15 also compares the variation in time with point P1 (the slope toe) displacement and the FOS of slope. In the early stage of rainfall (<2 h), the FOS only slightly decreases corresponding to insufficient accumulated pore-water pressure of slope. Between 2 h and 7 h, the FOS decreased monotonically, indicating that rising pore-water pressure and dissipating matric suction both decreased effective stresses simultaneously, and it corresponded to the rapid increase of displacement observed in P1, meaning that significant deformation was preceded by hydro-mechanical coupling. After about 7 h, the seepage process evolves to a quasi-steady state, in which the pore-water pressure and gravitation reach a transient equilibrium, and as such, the rate of displacement growth and the rate of decrease in FOS both become significantly smaller; the slope falls into the region of trivial or ultimate stability.

Figure 15

Figure 15. Evolution of monitoring point P1 displacement and slope FOS under non-uniform rainfall.

Figure 16 shows the changes of pore-water pressure at two monitoring points: P4 (in Rain-shielded area) and P5 (in Rain-exposed area). In initial period (0–2 h), the pore water pressure at P5 increases very fast and the matric suction sharply drops off, with an extremely large growth rate, and the pore water pressure of P4 increases slowly and the high suction remains. The difference between rainfall flux of the two zones creates a local high saturation point in the infiltration zone. This gives rise to a significant lateral hydraulic gradient from the infiltration zone to the adjacent shaded zone, favoring lateral infiltration of the moisture in infiltration zone to the shaded zones surrounding it. Over the 2–6 h period, the pore pressure at point P5 rises and tends to stabilize, but the speed of increasing pore pressure gets lower, while the speed of increasing pore pressure at P4 gets higher dramatically. The slopes of the lines at point P5 and P4 fluctuate between them, which means that after a lag time, the shaded area starts to see the effects of lateral seepage greatly reducing negative pore pressure. Unbalanced pore-pressure distributions and the lateral water flow interaction caused uneven suction dissipation in the middle slope region.

Figure 16

Figure 16. Pore-water pressure evolution at monitoring points P4 and P5 under non-uniform rainfall conditions.

As rainfall continues, the pore pressures at both two observation points are gradually reaching the steady state and the speed of their change gradually decreases, indicating that the shallow layer of the slope is gradually achieving a stable moisture state in the redistribution. At this time, the difference in pore pressures of point P4 and point P5 is only about 43.83 Pa, meaning that both points are close to having reached the same head. Overall, the spatial infiltration contrasts created by the non-uniform rain mainly influence the temporal variation of pore-pressure and the graduated hydromechanical response of the slope. In the early stage (0–6 h), the saturation evolution in the infiltration zone accelerates, and the shaded zone lags further, forming a distinguished ‘time lag difference’ in the two areas. However, at the later stage (6–24 h), redistribution effects decaying this difference become weaker, and finally, both regions reach a common hydraulic head.

Relative to uniform rainfall, non-uniform boundary conditions leave the final steady-state pore-pressure values essentially unchanged. In contrast, they exert a strong influence on the rate and timing of the pore-pressure increase. This early spatial difference causes a rapid drop of matric suction immediately below an infiltration area, with local reduction of effective stress and a reduction in shear strength. As lateral moisture diffusion proceeds, then other shaded zones also become affected, causing the total strength of the slope to drop in turn with a delay. Thus, the slope response when there is an uneven rain distribution is rather stage dependent and spatially concentrated, with more local effects on stability than in case of uniform rainfall. Uneven evolution of pore-water pressure has far more local effects on the local stability of the slope than uniform rainfall conditions.

Figure 17 presents the distribution of plastic strain and associated FOS at different times in the rainfall process. Before the initiation of rainfall (0–2 h), the slope is in a stable state. The FOS gradually decreases from 2.52 to 2.42. The larger infiltration flux in the infiltration zone quickly saturated the upper silty clay layer, so the effective stress locally reduced. Meanwhile, seepage convergence at the slope toe induces slight strain concentration. The deformations at this stage are still confined to the shallow layers and represent the early localized deformation induced by hydraulics. As rainfall continues (2–7 h), the pore pressure steadily accumulates, and lateral spread of pressure takes place along the interlayer weak zones, strengthening the downward diffusion of suction and further reduction of local strength. With the shear band tilting gradually downward along the interlayer weak zones and then turning into a continuous band like distribution, the largest shear strain rate changes sharply, and the peak deformation zone decreases towards the toe of the slope. This means that the slope has entered the stage from localized softening deformation into banded extension deformation transition. Interaction between strain concentration and pore pressure redistribution in this stage also results in coupled effects of infiltration and redistribution under non-uniform rainfall being the cause of localized extension. By 24 h, the shear band along the weak surface had become almost continuous and almost stable. Deformation was concentrated in the middle and lower slope and near the slope toe. Deformation occurred in the middle and lower slope and near the slope toe, more obvious shear strain rate peaks could be found; this deformation band outside the band area has a low strain rate. At this step, the value of FOS was reduced to 1.97, which indicates that the total shear strength of the slope largely decreases. The geometry of the shear band matches that of the interlayer weak zone and bends to upward at the slope toe, which indicates a shallower occurrence of translational failure than the simple sliding of the soil block and accompanied by internal tension break.

Figure 17

Figure 17. Evolution of maximum shear strain rate in the slope under non-uniform rainfall conditions.

Compared with the uniform rainfall case, the non-uniform scenario consistently shows a more concentrated pore-pressure increase beneath the shaded zone, a deeper wetting front in the convergence area, and an earlier decrease in FOS, indicating stronger localized instability.

6 Discussions

This study characterized rainfall infiltration–stability responses of slopes for non-uniform rainfall conditions due to complex surface features under mountain PV fields. The main goal of our work is to address the deficiency of ordinary uniform rainfall rain simulation models on the shading and runoff convergence effects of PV installations.

It was observed in rainfall redistribution (with PV panels) that serious spatiotemporal differences occur in infiltration intensity, leading to different seepage regions of injected water and different changes of pore-pressure evolution along the slope. The results show that this effect is not only affecting the propagation rate of the surface wetting front but is also facilitating local pore pressure buildup and unsynchronized reduction of slope stability. This localized pore-pressure intensification is consistent with the sensitivity analyses on unsaturated hydraulic behavior reported by Rahimi et al. (2011), who demonstrated that spatial variability in hydraulic properties or infiltration input can strongly amplify localized instability tendencies. Consistent with the local FOS field theory of Lu et al. (2012b), the slope failure occurring for arbitrary boundary loads proposed in the present study does not happen at the same time throughout the slope, but has the spatially differentiated features of triggering; that is, it does not trigger uniform spreading but rather the local instability prefers to develop in some parts of the slope. These findings indicate that neglecting rainfall heterogeneity in PV-covered slopes may underestimate local failure risk and misidentify the most vulnerable regions, particularly around infiltration convergence zones and shaded–unshaded transitions.

Similarly, uniform and non-uniform rain, the two series of FOS both present a similar three-stage evolution tendency in FOS over time, rapid deterioration, slow decay and accelerated decline (Figure 18). The three-stage FOS evolution observed in both rainfall patterns indicates that suction dissipation, transient hydraulic adjustment, and shear-band coalescence act sequentially as the dominant failure-controlling mechanisms. When identical rainfall duration is considered, the FOS value of non-uniform rainfall is ever slight (0.02–0.05) larger than that of the uniform rainfall. This is caused by limited infiltration under shaded areas and focused runoff within infiltration areas, leading to a temporal differentiating triggering rather than simultaneous degradation over the whole slope.

Figure 18

Figure 18. The temporal evolution of the FOS under uniform and non-uniform rainfall conditions.

The above approaches are simplified and should be improved in forthcoming studies. In the current study, the numerical model has a simplified three-layer stratigraphic structure and uses a fixed rainstorm intensity and duration and does not consider the influence of natural surface roughness and longitudinal runoff cooperation that both have a great impact on the local seepage condition.

Several limitations of this study should be acknowledged. The numerical model employs a simplified three-layer stratigraphy, fixed rainfall intensity, and idealized surface conditions, without explicitly representing micro-topography, longitudinal runoff cooperation, or crack-guided preferential flow, all of which may modify near-surface hydraulic connectivity and infiltration pathways. In addition, temperature–evaporation–moisture redistribution processes were not included, potentially affecting long-term pore-pressure evolution and leading to conservative estimates of suction recovery or slightly higher simulated FOS values compared with real conditions. The model also lacks field-monitored rainfall–pore-pressure–deformation data for direct validation, although engineering parameters, rainfall records, and PV panel geometry were taken from the actual project. Future work should integrate field-scale rainfall–infiltration experiments, long-term monitoring data, and parameter back-analysis to enhance model calibration and better capture the coupled hydro-mechanical response of slopes subjected to PV-induced non-uniform rainfall.

7 Conclusion

This study provides the first quantitative evaluation of the hydro-mechanical effects of PV-induced non-uniform rainfall on multi-layered slopes using a coupled FLAC3D framework, addressing a gap not covered in previous rainfall-induced landslide studies. To achieve this, this study developed a coupling hydro-mechanical modelling technique to study the rainfall-induced slope instability of mountain PV slope with spatially non-uniform infiltrations. Integrating Richards equation, Van Genuchten model and SRM all into a developed code in FLAC3D platform, we could well solve the variation of pore-water pressure, matric suction and shear strength over different (uniform and non-uniform) rainfalls.

The results demonstrate that PV panel shading, and surface runoff redistribution produce great temporal change in rainfall infiltration, causing a slowed and different relaxation of suction and uneven pore-pressure development in the slope. Compared with uniform rainfall, non-uniform rainfall led to localized saturation under the inflection region earlier, but the wetting below the shaded side later, leading to lateral hydraulic gradients and slight concentration of pore-pressure in a relatively narrow region near the toe of the slope.

The simulation results indicate that slope failure is governed by the coupling between non-uniform infiltration and interlayer hydraulic contrast. The stability evolves through three stages: (Ⅰ) a rapid-decline stage driven by the loss of matric suction; (Ⅱ) a transitional equilibrium stage controlled by transient hydraulic balance; and (Ⅲ) a final acceleration stage characterized by shear-band coalescence and failure.

The comparison study based on the numerical model proved that neglecting rainfall heterogeneity leads to an underestimation of localized failure risk in mountainous photovoltaic construction. Our modelling framework provides a validated computational tool for predicting slope responses under rainfall redistribution and partial shading induced hydrologic contrasts.

Data availability statement

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

Author contributions

DW: Writing – review and editing, Investigation. LT: Data curation, Writing – review and editing. QL: Writing – review and editing, Investigation. XY: Writing – review and editing. HH: Writing – review and editing. DT: Writing – review and editing. JY: Writing – review and editing. JZ: Writing – original draft, Formal Analysis, Conceptualization, Methodology, Writing – review and editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This work has been supported by the National Natural Science Foundation of China (NSFC) (No. 12402088), PetroChina Southwest Oil and Gas Field Company, (NO. 24H1697) and Sichuan University Interdisciplinary Innovation Fund (No. 0082204153389).

Acknowledgements

We would like to thank all individuals and organizations who contributed to this work, including technical and administrative support.

Conflict of interest

Authors DW, QL, LT, XY, HH and DT were employed by PetroChina Southwest Oil and Gas Field Company. Author JY was employed by Sichuan Huasheng Energy Development Group Co., Ltd.

The remaining author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Basher, L. R., and Ross, C. W. (2001). Role of wheel tracks in runoff generation and erosion under vegetable production on a clay loam soil at Pukekohe, New Zealand. Soil Tillage Res. 62 (3-4), 117–130. doi:10.1016/s0167-1987(01)00220-3

CrossRef Full Text | Google Scholar

Borja, R. I., and White, J. A. (2010). Continuum deformation and stability analyses of a steep hillside under rainfall infiltration. Acta Geotech. 5 (1), 1–14. doi:10.1007/s11440-009-0108-1

CrossRef Full Text | Google Scholar

Bryant, W. R. (2003). Permeability of clays, silty-clays and clayey-silts.

Google Scholar

Cai, F., and Ugai, K. (2004). Numerical analysis of rainfall effects on slope stability. Soils Found. 44 (2), 79–94. doi:10.1061/(asce)1532-3641(2004)4:2(69)

CrossRef Full Text | Google Scholar

Cao, L., Lv, M., Li, C., Sun, Q., Wu, M., Xu, C., et al. (2025). Effects of crosslinking agents and reservoir conditions on the propagation of fractures in coal reservoirs during hydraulic fracturing. Reserv. Sci. 1 (1), 36–51. doi:10.62762/rs.2025.494074

CrossRef Full Text | Google Scholar

Cho, S. E., and Lee, S. R. (2001). Instability of unsaturated soil slopes due to infiltration. Comput. Geotechnics 28 (3), 185–208. doi:10.1016/s0266-352x(00)00027-6

CrossRef Full Text | Google Scholar

Crosta, G. B., and Frattini, P. (2008). Rainfall thresholds for triggering soil slips and debris flow. Nat. Hazards Earth Syst. Sci. 8 (1), 31–50.

Google Scholar

Dai, F. C., and Lee, C. F. (2002). Landslide characteristics and slope instability analysis along the Lantau north coastal area, Hong Kong. Eng. Geol. 63 (3-4), 251–267.

Google Scholar

Fan, L., Lehmann, P., and Or, D. (2015). Effects of hydromechanical loading history and antecedent soil mechanical damage on shallow landslide triggering [J]. Journal of Geophysical Research: Earth Surface 120 (10), 1990–2015.

CrossRef Full Text | Google Scholar

Fleckenstein, J. H., Krause, S., Hannah, D. M., and Boano, F. (2010). Groundwater-surface water interactions: new methods and models to improve understanding of processes and dynamics. Adv. Water Resour. 33, 1291–1295. doi:10.1016/j.advwatres.2010.09.011

CrossRef Full Text | Google Scholar

Fredlund, D. G., and Rahardjo, H. (1993a). Soil mechanics for unsaturated soils. New York: Wiley.

Google Scholar

Fredlund, D. G., and Rahardjo, H. (1993b). Soil mechanics for unsaturated soils. John Wiley and Sons.

Google Scholar

Garcia Aristizabal, E. F., and Riveros Jerez, C. A. (2011). Builes Brand M A. Influence of rainfall intensity on infiltration and deformation of unsaturated soil slopes [J]. Dyna 78 (170), 116–124.

Google Scholar

Garcia-Bengochea, I., Altschaeffl, A. G., and Lovell, C. W. (1979). Pore distribution and permeability of silty clays. J. Geotechnical Engineering Division 105 (7), 839–856. doi:10.1061/ajgeb6.0000833

CrossRef Full Text | Google Scholar

Iverson, R. M. (2000). Landslide triggering by rain infiltration. Water Resour. Res. 36 (7), 1897–1910. doi:10.1029/2000wr900090

CrossRef Full Text | Google Scholar

Jing, J., Hou, J., Sun, W., Chen, G., Ma, Y., and Ji, G. (2022). Study on influencing factors of unsaturated loess slope stability under dry-wet cycle conditions. J. Hydrology 612, 128187. doi:10.1016/j.jhydrol.2022.128187

CrossRef Full Text | Google Scholar

Jordán, A., and Martínez-Zavala, L. (2008). Soil loss and runoff rates on unpaved forest roads in southern Spain after simulated rainfall. For. Ecol. Manag. 255 (3-4), 913–919. doi:10.1016/j.foreco.2007.10.002

CrossRef Full Text | Google Scholar

Kalantari, Z., Briel, A., Lyon, S. W., Olofsson, B., and Folkeson, L. (2014). On the utilization of hydrological modelling for road drainage design under climate and land use change. Sci. Total Environ. 475, 97–103. doi:10.1016/j.scitotenv.2013.12.114

PubMed Abstract | CrossRef Full Text | Google Scholar

Kristo, E., Rahardjo, H., and Leong, E. C. (2017). Effect of variations in rainfall intensity on slope stability in Singapore. Eng. Geol. 228, 234–244. doi:10.1016/j.iswcr.2017.07.001

CrossRef Full Text | Google Scholar

Li, L., Ma, S., Liu, X., Liu, J., Lu, Y., Zhao, P., et al. (2024). Coal measure gas resources matter in China: review, challenges, and perspective. Phys. Fluids 36 (7), 071301. doi:10.1063/5.0218328

CrossRef Full Text | Google Scholar

Li, M., Liu, J., and Xia, Y. (2025). Risk prediction of gas hydrate formation in the wellbore and subsea gathering system of deep-water turbidite reservoirs: case analysis from the South China Sea. Reserv. Sci. 1 (1), 52–72. doi:10.62762/rs.2025.567907

CrossRef Full Text | Google Scholar

Li, L., Abdallah, K. B., Hamdi, E., Hou, B., Cui, Z., Elsworth, D., et al. (2026). CO2-Enhanced multiphase flow in heterogenous coal measures: thermal-hydraulic-mechanic (THM) model for enhancing gas Co-Production with CO2 geo-sequestration. Fuel 407, 137479. doi:10.1016/j.fuel.2025.137479

CrossRef Full Text | Google Scholar

Liu, Y., Zhao, Y., Vanapalli, S. K., and Mehmood, M. (2024). Soil-water characteristic curve of expansive soils considering cumulative damage effects of wetting and drying cycles. Eng. Geol. 339, 107642. doi:10.1016/j.enggeo.2024.107642

CrossRef Full Text | Google Scholar

Lu, N., and Likos, W. J. (2004). Unsaturated soil mechanics. New York: Wiley.

Google Scholar

Lu, N., Godt, J. W., and Wu, D. T. (2012a). Analysis of rainfall-induced slope instability using a field of local factor of safety. Water Resources Research, 48, W09524. Rahimi A, Rahardjo H, Leong E C. Effect of hydraulic properties of soil on rainfall-induced slope failure. Eng. Geol. 114 (3-4), 135–143. doi:10.1016/j.enggeo.2010.04.010

CrossRef Full Text | Google Scholar

Lu, N., Şener-Kaya, B., Wayllace, A., and Godt, J. W. (2012b). Analysis of rainfall-induced slope instability using a field of local factor of safety. Water Resour. Res. 48 (9). doi:10.1029/2012wr011830

CrossRef Full Text | Google Scholar

Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12, 513–522. doi:10.1029/wr012i003p00513

CrossRef Full Text | Google Scholar

Ng, C. W. W., and Shi, Q. (1998). A numerical investigation of the stability of unsaturated soil slopes subjected to transient infiltration. Comput. Geotechnics 22 (1), 1–19. doi:10.1016/S0266-352X(97)00036-0

CrossRef Full Text | Google Scholar

Nguyen, T. S., Likitlersuang, S., Tanapalungkorn, W., Phan, T. N., and Keawsawasvong, S. (2022). Influence of copula approaches on reliability analysis of slope stability using random adaptive finite element limit analysis. Int. J. Numer. Anal. Methods Geomechanics 46 (12), 2211–2232. doi:10.1002/nag.3385

CrossRef Full Text | Google Scholar

Nie, W., Li, C., Hu, J., Saffari, P., Wang, W., and Luo, M. (2023). Spatial variation of physical and mechanical properties of tailings under different rainfall intensities and the interaction pattern. Geomechanics Geophys. Geo-Energy Geo-Resources 9 (1), 86. doi:10.1007/s40948-023-00625-0

CrossRef Full Text | Google Scholar

Oh, S., and Lu, N. (2015). Slope stability analysis under unsaturated conditions: case studies of rainfall-induced failure. Can. Geotechnical J. 52 (12), 1–14. doi:10.1016/j.enggeo.2014.11.007

CrossRef Full Text | Google Scholar

Rahimi, A., Rahardjo, H., and Leong, E. C. (2010). Effect of hydraulic properties of soil on rainfall-induced slope failure. Comput. Geotechnics 37 (3), 313–320. doi:10.1016/j.enggeo.2010.04.010

CrossRef Full Text | Google Scholar

Rahimi, A., Rahardjo, H., and Leong, E. C. (2011). Effect of antecedent rainfall patterns on rainfall-induced slope failure. Eng. Geol. 120 (1-4), 1–14. doi:10.1061/(asce)gt.1943-5606.0000451

CrossRef Full Text | Google Scholar

Richards, L. A. (1931). Capillary conduction of liquids through porous mediums. Physics 1, 318–333. doi:10.1063/1.1745010

CrossRef Full Text | Google Scholar

Rosly, M. H., Mohamad, H. M., Bolong, N., and Harith, N. S. H. (2023). Relationship of rainfall intensity with slope stability. Civ. Eng. J. 9 (Special Issue), 75–83. doi:10.28991/cej-sp2023-09-06

CrossRef Full Text | Google Scholar

Song, Y. S., Chae, B. G., Kim, K. S., Park, J. Y., Oh, H. J., and Jeong, S. W. (2021). A landslide monitoring system for natural terrain in Korea: development and application in hazard evaluations. Sensors 21 (9), 3040. doi:10.3390/s21093040

PubMed Abstract | CrossRef Full Text | Google Scholar

Sophocleous, M. (2002). Interactions between groundwater and surface water: the state of the science. Hydrogeol. 10, 52–67. doi:10.1007/s10040-001-0170-8

CrossRef Full Text | Google Scholar

Subedi, R., Aryal, M., Dawadi, A., and Bastola, B. (2025). Seepage and slope stability analysis of a reservoir side slope: a numerical modeling approach.

Google Scholar

Terzaghi, K., Peck, R. B., and Mesri, G. (1996). Soil mechanics in engineering practice [M]. John wiley & sons 68 (142), 149–150.

Google Scholar

Tu, X. B., Kwong, A. K. L., Dai, F. C., Tham, L., and Min, H. (2009). Field monitoring of rainfall infiltration in a loess slope and analysis of failure mechanism. Eng. Geol. 105 (3-4), 134–150. doi:10.1016/j.enggeo.2008.11.011

CrossRef Full Text | Google Scholar

Van Genuchten, M. T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898. doi:10.2136/sssaj1980.03615995004400050002x

CrossRef Full Text | Google Scholar

Wang, J. H., Xu, W. J., and Liu, X. X. (2024). A slope stability analysis method considering the rainfall hydrology process. Eng. Geol. 343, 107775. doi:10.1016/j.enggeo.2024.107775

CrossRef Full Text | Google Scholar

Wang, K., Zhang, L., and Yu, M. (2025). Three-dimensional coupled hydrological and strength-reduction analysis for rainfall-induced slope failure prediction. Front. Earth Sci. 13, 1522876.

Google Scholar

Wu, L. Z., Huang, R. Q., Xu, Q., Zhang, L. M., and Li, H. L. (2015a). Analysis of physical testing of rainfall-induced soil slope failures. Environ. Earth Sci. 73, 8519–8531. doi:10.1007/s12665-014-4009-8

CrossRef Full Text | Google Scholar

Wu, L. Z., Huang, R. Q., and Xu, Q. (2015b). Experimental study on rainfall infiltration and slope failure mechanism. Environ. Earth Sci. 73, 8520–8525.

Google Scholar

Yuan, C., Li, Q., Nie, W., and Ye, C. (2023). A depth information-based method to enhance rainfall-induced landslide deformation area identification. Measurement 219, 113288. doi:10.1016/j.measurement.2023.113288

CrossRef Full Text | Google Scholar

Zhang, L. L., Zhang, L. M., and Tang, W. H. (2005). Rainfall-induced slope failure considering variability of soil properties [J]. Geotechnique 55 (2), 183–188.

CrossRef Full Text | Google Scholar

Zhang, L. M., Zhou, Y., and Gao, L. (2014). Rainfall-induced slope failures in China: mechanisms and mitigation. Eng. Geol. 179, 97–111.

Google Scholar

Zhang, S. Q., Abdallah, K. B., Li, L., Hamdi, E., and Liu, J. (2025a). Multiphase flow controlled by synergistic injection-production pressure: enabling CO2 geo-sequestration with additional gas recovery from vertically heterogeneous depleted shale reservoirs. Fuel 398, 135591. doi:10.1016/j.fuel.2025.135591

CrossRef Full Text | Google Scholar

Zhang, M., Wang, W., Hao, Y., Liu, X., Zhang, R., Tian, Z., et al. (2025b). Study on the formation mechanisms and stability of a gently inclined and shallow accumulative landslide. Sci. Rep. 15 (1), 15246. doi:10.1038/s41598-025-98031-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, C., and Sun, Q. (2019). Global rainfall-induced landslide database and analysis. Landslides 16 (1), 1–14.

Google Scholar

Zhou, T., Wang, J., and Xu, D. (2024). Quantitative evaluation of rainfall pattern and soil hydromechanical coupling in slope failure thresholds. Front. Geohazards Georisks 3, 1432057.

Google Scholar

Zienkiewicz, O. C., Chan, A. H. C., Pastor, M., Paul, D. K, and Shiomi, T. (1990). Static and dynamic behaviour of soils: a rational approach to quantitative solutions. I. Fully saturated problems. Proc. R. Soc. Lond. A. Math. Phys. Sci. 429 (1877), 285–309.

CrossRef Full Text | Google Scholar