A Hybrid Prediction Model for Rock Reservoir Bank Slope Deformation Considering Fractured Rock Mass Parameters

Abstract

Deformation monitoring data provide a direct representation of the structural behavior of reservoir bank rock slopes, and accurate deformation prediction is pivotal for slope safety monitoring and disaster warning. Among various deformation prediction models, hybrid models that integrate field monitoring data and numerical simulations stand out due to their well-defined physical and mechanical concepts, and their ability to make effective predictions with limited monitoring data. The predictive accuracy of hybrid models is closely tied to the precise determination of rock mass mechanical parameters in structural numerical simulations. However, rock masses in rock slopes are characterized by intersecting geological structural planes, resulting in reduced strength and the creation of multiple fracture flow channels. These factors contribute to the heterogeneous, anisotropic, and size-dependent properties of the macroscopic deformation parameters of the rock mass, influenced by the coupling of seepage and stress. To improve the predictive accuracy of the hybrid model, this study introduces the theory of equivalent continuous media. It proposes a method for determining the equivalent deformation parameters of fractured rock mass considering the coupling of seepage and stress. This method, based on a discrete fracture network (DFN) model, is integrated into the hybrid prediction model for rock slope deformation. Engineering case studies demonstrate that this approach achieves a high level of prediction accuracy and holds significant practical value.

1. Introduction

Reservoirs play a crucial role in regulating water supply, flood control, disaster reduction, and ecological conservation. However, the safe operation of dams and reservoir bank slopes faces significant challenges due to various adverse factors such as water pressure, heavy rainfall, and extreme weather events [1]. To prevent engineering accidents, it is essential to monitor and diagnose the safety status of reservoir bank slopes [2,3,4,5,6]. Monitoring data provide the most comprehensive and intuitive reflection of the evolutionary process of slope structural behavior under multiple factors. Currently, two primary methodologies are employed for monitoring, data analysis, and prediction.

The first approach involves numerical simulation methods such as finite element models and geological mechanics models to simulate the behavior characteristics of slopes under various influencing factors, predicting slope deformation and stability. Zhuang et al. [7] investigated the impact of soft interlayers, faults, and extreme rainfall on the stability of high rock slopes, combining numerical analysis with strength reduction methods to calculate landslides’ stability under natural and extreme rainfall conditions, validated with field monitoring data. Paswan [8] developed a rainfall-induced slope monitoring system based on Micro-Electro–Mechanical System (MEMS) sensors. demonstrating effective monitoring of both gradual and sudden slope movements induced by rainfall. However, due to the highly nonlinear nature of geological conditions, the complex properties of rock masses involved in predicting rock slope deformation and instability timing, and the intricacies of mechanics [9], the consideration of rheological or geometric nonlinear numerical simulation methods currently struggles to comprehensively reflect the dynamic evolution of rock slope deformation under various complex conditions.

The second approach utilizes prediction models based on monitoring data and various mathematical theories, such as statistical analysis, fuzzy mathematics, and artificial intelligence theories. For instance, Tonini [10] first introduced the idea that factors affecting dam displacement included water pressure components, temperature components, and time-dependent components, representing water pressure and temperature components using third-degree polynomials. To address the issue of noise contamination in prototype data of dam foundation seepage pressure, Zhu et al. [11] employed the Variational Mode Decomposition (VMD)–Wavelet Packet Denoising method to denoise the noisy prototype data, develo** an improved Time-Series Convolutional Network (ITCN) model for effective prediction of dam foundation seepage pressure data.

Lin et al. [12,13] decomposed deformation sequences based on temporal characteristics and employed various artificial intelligence algorithms to predict each sequence component, yielding favorable prediction results. Dong [14] introduced a real-time wireless monitoring system to obtain high-frequency overall data describing the deformation characteristics of steep slopes, using a Deep Autoregressive Network (DeepAR) model to predict slope displacements and assess the safety of excavated slopes. Du et al. [15] integrated weaker learners, proposing a novel ensemble learner to mine field data acquired from Ground-Based Synthetic Width Radar (GB-SAR) to create a slope deformation prediction model. However, mathematical models are primarily suited for datasets with longer sequences of data and sometimes lack the ability to provide a fundamental mechanical understanding of slope deformation behavior.

Addressing these issues, Bonaldi et al. [16] combined the strengths of empirical data and numerical simulations, proposing a hybrid model for the deformation of gravity dams. In this model, primary deformation components or some deformation components are calculated using finite element numerical analysis, referred to as deterministic components with clear physical concepts. The remaining components are fitted based on statistical models using monitoring data. This approach proves effective for monitoring and prediction, especially when monitoring data are limited. However, the accuracy of the hybrid model’s predictions depends heavily on the selection of rock mass mechanical parameters during the calculation of deterministic components. Rock reservoir bank slopes are often situated in mountainous and complex geological settings, with internal rock masses segmented by various geological structural planes, resulting in lower strengths and the formation of fracture permeation pathways. These factors result in the performance parameters of rock masses, such as their macroscopic deformation, displaying non-uniformity, anisotropy, and size effects. Additionally, numerous fractures within rock reservoir bank slopes constitute primary seepage channels, with continuous external water environments altering fracture widths and deformation responses. Studying the mechanical parameters of rock reservoir bank slopes requires considering the size effects of these parameters and the coupled effects of the seepage field and stress field.

In conclusion, this paper considers the influence of geological conditions and seepage factors on the mechanical parameters, deformation behavior, and monitoring data of slopes. It proposes a method for determining the equivalent deformation parameters of fractured rock masses under seepage–stress coupling. Subsequently, the deformation behavior of slopes under seepage conditions is computed based on the theory of equivalent continuous media, and a deterministic model is established to determine rainfall components. Finally, a hybrid model for slope deformation prediction is constructed by incorporating temperature components and time-dependent factors. The research roadmap is illustrated in Figure 1.

2. Hybrid Prediction Model for Rock Reservoir Bank Slope Deformation

Rainfall, humidity, evaporation, temperature, groundwater, and snowmelt all influence slope deformation [17]. These influencing factors often number in the tens or even dozens. Including all these factors in a regression equation can lead to the degradation of the coefficient matrix, making it unsolvable or reducing the precision. Therefore, before establishing the regression equation, researchers often use methods such as the Maximal Information Coefficient (MIC) method [18], stepwise regression [19], and sensitivity analysis [20] to select the factors that significantly contribute to slope deformation. This ensures that the resulting regression equation includes only significant factors and excludes insignificant ones.

This paper introduces a multifactor regression model based on the stepwise regression method. At each step, F-tests and Variance Inflation Factor (VIF) tests are performed to add or remove factors, ensuring that only significant factors are included in the regression equation and reducing multicollinearity among the model variables. Through the analysis of deformation monitoring data from several reservoir slopes and feedback from actual engineering projects in China, the model identifies rainfall, temperature, and aging as the main factors influencing slope deformation $\delta$ (including cumulative horizontal and vertical displacements) [21,22]. Building on this model, the coefficients for the rainfall component ${\widehat{\delta }}_U$ are determined by fitting the equation to slope deformation data obtained from multiple numerical calculations. The coefficients for the temperature component ${\delta }_T$ and time-dependent component ${\delta }_{\theta }$ are determined by fitting the equation to observed slope deformation data. This results in the formulation of a hybrid model for predicting rock slope deformation, as shown in Equation (1).

$$\delta ={\widehat{\delta }}_U+{\delta }_T+{\delta }_{\theta }$$

(1)

① Rainfall Component ${\widehat{\delta }}_U$

The rainfall component is represented using a deterministic expression obtained through finite element analysis. Rainfall impacts slope deformation by infiltrating the ground, raising the groundwater level, and altering the slope’s water content. This process exhibits a time lag effect between slope deformation and rainfall [23,24,25]. Previous studies have shown that rainfall occurring more than 15 days prior can be neglected to avoid redundant calculations due to minimal impact [22]. Therefore, the observed rainfall within the 15 days prior to the observation date is used as the boundary condition for the finite element analysis in this paper. By fitting the deformation values of measurement points obtained from finite element analysis under various rainfall conditions into the deterministic model expression for the rainfall component, the rainfall fitting coefficients can be determined. There are several methods to define the specific expressions for rainfall components [26]. This paper adopts the approach from reference [27], selecting multiple antecedent average rainfall amounts as the rainfall factors:

$${\widehat{\delta }}_U={\displaystyle \sum _{i=1}^5\left(a_i{U}_i\right)}$$

(2)

where $U_i$ is the mean antecedent observed rainfall on the observed days, where $U_1, U_2, U_3, U_4$ is the mean rainfall of the antecedent 1, 2, 3, and 7 days, respectively. $a_i$ is the regression coefficient of the rainfall factor. The main steps for determining coefficients $a_i$ are as follows:

Step 1: Based on long-term rainfall observation data for the region where the slope is located, determine the range of daily rainfall to be used as boundary conditions;

Step 2: Randomly combine the daily rainfall values within this range to generate multiple 15-day rainfall sample sets;

Step 3: Using the input sample sets, perform numerical simulations to calculate the deformation at each measurement point on the slope. This results in training datasets composed of the influencing factor samples and the corresponding deformations at each measurement point;

Step 4: Using the obtained training dataset, fit the coefficients $a_i$ in Equation (2).

② Temperature Component ${\delta }_T$

Temperature variations can affect the fracture width and stresses in rock slopes, which consequently influence the stability of slopes. Unlike structures such as dams, where thermometers can be embedded during construction, temperature sensors are rarely installed within slopes. In cases where monitoring data are insufficient, temperature components are typically simulated using the periodic terms of sine and cosine functions. The formula is as follows [22,28]:

$${\delta }_T={\displaystyle \sum _{i=1}^m[b_{1i}(\sin \frac{2\pi it}{365}-\sin \frac{2\pi i{t}_0}{365})+b_{2i}(\cos \frac{2\pi it}{365}-\cos \frac{2\pi i{t}_0}{365})]}$$

(3)

where $t$ is the cumulative number of days from the observation date to the starting date, $t_0$ is the cumulative number of days from the first monitoring date of the modeling dataset to the initial monitoring date, and $b_{1i}$ and $b_{2i}$ are statistical coefficients for temperature factors, such as yearly ($i=1$) or semi-yearly ($i=2$) cycles. The coefficients $b_{1i}$ and $b_{2i}$ for the rainfall component, along with the coefficients $c_1$ and $c_2$ for the time-dependent component, are determined using regression equations. This is accomplished by subtracting the rainfall component from the observed deformation values at each measurement point and analyzing the relationship between the remaining deformation ${\delta }_{T+\theta }=\delta -{\widehat{\delta }}_U$ and influencing factors.

③ Time-Dependent Component ${\delta }_{\theta }$

The time-dependent component reflects the trend changes in slope deformation over time. It includes the deformation caused by the deterioration and fragmentation of the rock mass under the long-term action of water–rock, as well as the creep and plastic deformation of the rock mass and rock discontinuity. Considering the nonlinear characteristics of the time-dependent component, a combination of linear and nonlinear terms is used for slopes that have not yet entered the accelerated deformation stage [22,29]:

$${\delta }_{\theta }=c_1(\theta -{\theta }_0)+c_2(\ln \theta -\ln {\theta }_0)$$

(4)

where $\theta =t/100$, $t$ denotes the cumulative number of days from the observation date to the starting date, and $c_1$ and $c_2$ are regression coefficients for time-dependent factors. $c_1(\theta -{\theta }_0)$ corresponds to the linear portion of the rock creep curve, and $c_2(\ln \theta -\ln {\theta }_0)$ corresponds to the nonlinear portion. The method for determining the coefficients $c_1$ and $c_2$ for the time-dependent component is the same as that for determining the temperature component coefficients $b_{1i}$ and $b_{2i}$.

④ Expression of the Hybrid Model

In summary, considering the characteristics of rock reservoir bank slopes and accounting for the influence of initial measurements, the hybrid model for slope deformation monitoring data can be expressed as follows:

$$\begin{array}{l}\delta (t)=a_0+{\displaystyle \sum _{i=1}^n\left(a_i{U}_i\right)}+{\displaystyle \sum _{i=1}^m[b_{1i}(\sin \frac{2\pi it}{365}-\sin \frac{2\pi i{t}_0}{365})+b_{2i}(\cos \frac{2\pi it}{365}-\cos \frac{2\pi i{t}_0}{365})]}\\ +c_1(\theta -{\theta }_0)+c_2(\ln \theta -\ln {\theta }_0)\end{array}$$

(5)

where $a_0$ represents the constant term.

The hybrid model developed here encapsulates the combined effects of rainfall, temperature, time-dependent factors, and initial measurements to facilitate the prediction of deformations in rock reservoir bank slopes.

3. Determination of Rock Mass Mechanical Parameters

As shown in Figure 2, the mechanical properties of fractured rock masses vary with the size of the rock blocks. However, when the size reaches a certain value known as the mechanical Representative Elementary Volume (REV), the mechanical parameters of the rock mass become stable. The concept of the REV plays a pivotal role in the selection of mechanical models and parameter values across various phases of rock engineering, including surveying, testing, design, and construction. Initially introduced by Bear [30] in the context of groundwater flow in soils, the concept of REV has subsequently been applied to the study of rock masses. Shlomo [31] utilized in situ measurements to establish the REV of rock masses, while Kulatilake [32] employed finite element analysis to assess the representation of rock mass units. Ki Bok et al. [33] conducted an analysis of the REV of jointed rock mass using the Discrete Element Method (DEM). It is important to note that this study focuses only on rock mass with non-persistent joints. Persistent joints, which could influence the whole slope deformation, should be analyzed separately.

At present, various methods are utilized to determine the mechanical parameters of materials within rock reservoir bank slopes. Analytical methods typically assume regular, non-intersecting structural surfaces, making it challenging to accurately represent the actual structural characteristics of the rock mass. Experimental techniques provide parameters that only reflect the characteristics of rock masses within a limited range around rock samples or test points, failing to represent the macroscopic properties of the measured rock mass. Numerical methods, on the other hand, have gained prominence in the study of the REV of rock mass. They enable the reconstruction of discontinuity surfaces in rock masses by integrating geological survey data and stochastic sampling techniques. These methods facilitate the assignment of specific constitutive relationships or properties to discontinuity surfaces or rock blocks, allowing for numerical simulations of stress–strain relationships on these surfaces. In this study, finite element analysis is employed to determine the mechanical parameters of slope rock masses. By introducing complex joint elements and accounting for variations in joint widths, iterative calculations of the seepage field and stress field are conducted, enabling the determination of deformations in both rock media and fracture media within rock mass analysis units. Subsequently, equivalent deformation parameters for rock masses are computed based on the principle of deformation equivalence.

3.1. Numerical Simulation Analysis Method Considering Seepage–Stress Coupling

The finite element analysis software ABAQUS 6.14-4 is used to simulate the discrete fracture network (DFN). The governing equations and finite element solution models are thoroughly described in reference [32]. The block elements of the equivalent rock mass adopt the Mohr–Coulomb constitutive model, while the joint surfaces consist of composite elements, which include equal-thickness Goodman’s joint elements [34] and line elements. The relationship between mechanical joint widths and stress is represented using a parabolic equation [35].

3.1.1. Composite Fracture Elements

When conducting deformation analysis of fracture media considering seepage effects, two sets of node and element information come into play: equal-thickness Goodman’s joint elements for stress analysis and line elements for seepage analysis. As depicted in Figure 3, these two sets of information employ identical element numbers to denote the same computational domains, both corresponding to fractures, ensuring accurate transmission of all coupled parameters across all elements during the computational steps. It is worth noting that the line elements shown in Figure 3 do not possess thickness. However, for the computation of fracture permeability matrices, their thickness is assigned based on the joint widths at the current time step during the coupling calculation with the stress field.

3.1.2. Fracture Width Variation

Fracture width variation is fundamental in the coupled analysis of rock mass permeability and stress. Under specific normal stress conditions, changes in fracture flow are primarily attributed to changes in fracture width. The relationship between mechanical width variation and stress is established using existing empirical deformation formulas, facilitating the derivation of the link between equivalent hydraulic width and stress. This approach, grounded in extensive experimental research and practical applications, offers a robust theoretical basis and exposes the core impact of stress on fracture flow. Specifically, the relationship between mechanical width variation and stress is expressed using the following parabolic equation [35]:

$$\Delta V=\frac{{\sigma }_n}{k_{n0}+{\sigma }_n/b_{m0}}$$

(6)

where $\Delta V$ is mechanical width variation, $k_{n0}$ is the initial normal stiffness coefficient, ${\sigma }_n$ is the normal effective stress on the fracture surface, and $b_{m0}$ is the initial fracture width at zero normal stress. The above equation can be simplified as follows [35]:

$$\Delta V=\frac{{\sigma }_n}{k_n}$$

(7)

where $k_n=k_{n0}+{\sigma }_n/b_{m0}$ is the secant normal stiffness coefficient.

Considering that, under high-stress conditions, fracture surfaces typically do not fully close, thereby maintaining permeable pathways, a smaller residual fracture width should be assigned when the calculated mechanical deformation of the fracture surface equals or exceeds the initial fracture width. At this juncture, the relationship between fracture width and normal stress on the fracture is represented as follows [35]:

$$b_m=f\left({\sigma }_n\right)=\{\begin{array}{ll}{b}_{m0}\left(1-\frac{{\sigma }_n}{b_{m0}{k}_n}\right)\hfill & \left({\sigma }_n<b_{m0}{k}_n\right)\hfill \\ \mathrm{Res}\hfill & \left({\sigma }_n\ge b_{m0}{k}_n\right)\hfill \end{array}$$

(8)

where $b_m$ is the fracture width, and $\mathrm{Res}$ is the residual fracture width.

To apply Darcy’s law for fracture flow calculations, the mechanical fracture width derived from Equation (8) must be converted to an equivalent hydraulic fracture width, ensuring an accurate representation of stress’s impact on fracture flow. Additionally, since natural fracture surfaces typically possess roughness, a correction to the cubic law is necessary for single-fracture flow. Based on empirical data, Barton proposed a correlation between equivalent hydraulic width, mechanical fracture width, and the Joint Roughness Coefficient (JRC) [36]:

$$b_h=\frac{b_m^2}{JR{C}^{2.5}}$$

(9)

3.1.3. Coupled Analysis Process of DFN for Flow and Stress

By incorporating the obtained fracture width changes into the iterative computations of both the flow field and stress field, the final flow field and deformation field can be determined. The analysis process is illustrated in Figure 4.

3.2. Calculation of Equivalent Deformation Parameters Considering Seepage–Stress Coupling

As previously discussed regarding the Representative Elementary Volume (REV) concept, when the size of the equivalent continuous body exceeds the REV, the equivalent parameters that characterize its behavior stabilize. To obtain representative equivalent deformation parameters, finite element analysis models are established based on a stochastic fracture network, utilizing the fracture distribution characteristics of the slope rock mass obtained from on-site engineering surveys. Numerical simulation tests are conducted on finite element models of varying sizes and orientations to calculate the equivalent deformation parameters. The calculation process is illustrated in Figure 5.

4. Case Study

This section analyzes the equivalent mechanical parameters of fractured rock masses and the hybrid prediction model for the deformation of rock slopes at ZN Pumped Storage Power Station in China.

4.1. Overview of the Project

ZN Pumped Storage Power Station is located in Fujian Province, China, and boasts an installed capacity of 1200 MW, classifying it as a Type I large hydropower project. The infrastructure of the power station encompasses the upper reservoir, lower reservoir, water conveyance system, underground powerhouse, and other structures. The lower reservoir is nestled in a creek valley with roughly symmetrical dam surroundings. The area surrounding the reservoir exhibits steep terrain, characteristic of highly dissected, low-to-mid-elevation mountain slopes, forming a deep “V” shape in the valley. These mountainous terrains feature steep slopes ranging from 38° to 60°. The geological conditions at the dam site are intricate, marked by several secondary faults, rock veins, and unfavorable geological structures. Some natural slopes in the area exhibit deep unloading rock masses. The exposed bedrock formations in the lower reservoir (dam) area primarily consist of the Upper Nan-yuan Formation from the Late Jurassic period, late-stage Yanshan intrusive rocks, and Quaternary surface formations. The predominant rock type is potassic feldspar granite, composed of minerals like potassium feldspar, plagioclase, quartz, and biotite. This granite exhibits a medium-to-fine-grained granitic structure with blocky features. Both sides of the reservoir area feature bedrock outcrops with varying degrees of weathering. The full weathered layer’s thickness is limited, and various diabase dikes (βμ) and a quartz diorite porphyry dike (ηoπ) can be found.

Our study focuses on the rock slope segment stretching from 0 + 114.5 to 0 + 162.5 m downstream of ZN Pumped Storage Power Station’s lower reservoir, with elevations ranging from 190 to 280 m. Based on the information provided by the design department, the positions of the monitoring instruments and the geological cross-section of the typical profile are shown in Figure 6.

4.2. Establishment of Numerical Models

(1) Distribution characteristics of slope joints

Below, we focus specifically on the Mild Weathering (Upper Section) Potassium Feldspar Granite (Rock Mass Rating III2A) as an example, detailing the entire process from stochastic fracture network modeling to the calculation of final equivalent deformation parameters. The potassic feldspar granite in this area exhibits a blocky-to-sub-blocky pattern. Extensive joint development is observed, with some joints extending for longer distances and featuring flat surfaces. These fracture surfaces often display iron–manganese staining. Geological map** and tunnel data examination revealed that the major structural discontinuity groups on the left and right banks consist of three sets each, characterized by the major structural discontinuities with suggested combinations as presented in

Table 1. Slope structural surface combination characteristics table.

Location	Fracture Set	Strike	Characteristics	Spacing (cm)
Left Bank	1	${\mathrm{N}50~60}^{\circ }\mathrm{E} \mathrm{NW}\left(\mathrm{SE}\right)\angle 7{0~90}^{\circ }$	Long extensions, flat surfaces, multiple surfaces with iron–manganese staining	20~50
	2	${\mathrm{N}60~70}^{\circ }\mathrm{E} \mathrm{NWZ}\angle 8{0~90}^{\circ }$	Long extensions, flat surfaces, some iron staining or mud infill	30~50
	3	${\mathrm{N}40~50}^{\circ }\mathrm{E} \mathrm{NW}\left(\mathrm{SE}\right)\angle 3{0~40}^{\circ }$	Long extensions, flat surfaces, iron staining	20~50
Right Bank	1	${\mathrm{N}80~90}^{\circ }\mathrm{E} \mathrm{NWZ}\angle {80~90}^{\circ }$	Long extensions, flat surfaces, iron staining	20~50
	2	${\mathrm{N}70~80}^{\circ }\mathrm{E} \mathrm{NWZ}\angle {70~90}^{\circ }$	Long extensions, flat surfaces, iron staining	10~50
	3	${\mathrm{N}60~70}^{\circ }\mathrm{E} \mathrm{NWZ}\angle {60~90}^{\circ }$	Long extensions, flat surfaces, iron–manganese staining	20~50

.

(2) Construction of Finite Element Models for Slope Rock Mass

To determine the REV of the slope rock mass, we studied variations in equivalent parameters of rock mass analysis units of various sizes and orientations. Firstly, it is necessary to establish a stochastic fracture network. The modeling principle relies on the statistical parameters and patterns of on-site fractures, which are then simulated using computer visualization techniques. The modeling process primarily involves the following steps:

Step 1: Geological surveys of the excavation face are conducted to collect data on rock mass structural features;

Step 2: Based on the structural orientation characteristics, the data are categorized and subjected to statistical analysis for each category;

Step 3: Probability distribution models representing structural feature parameters are constructed. In these models, dip angles follow a normal distribution, while trace lengths and spacings follow negative exponential distributions;

Step 4: The Monte Carlo method is applied to perform random sampling of fracture feature parameters, generating a series of random variables conforming to the probability models mentioned above. These variables are then used to construct the rock mass fracture network model;

Step 5: The fracture parameters in the model are validated and adjusted through comparison with field survey data.

After utilizing the fracture distribution pattern described in Table 1 and employing the Monte Carlo method, we created a random fracture network. This network was then imported into the Hypermesh 2023 software. Various square rock mass analysis unit models were constructed with different sizes and orientations, using the center point of the fracture network model as the reference point, as shown in Figure 7.

According to the research outlined in reference [37], the maximum model size should be at least three times the average trace length, while the minimum size should be similar to the average trace length. As such, the side lengths of models in various directions were set as follows: 4 m, 6 m, 8 m, 10 m, 12 m, 14 m, 16 m, and 18 m. The directions for each model were set at 0°, 30°, 60°, 90°, 120°, and 150°. Additionally, rock mass analysis units in directions 180°, 210°, 270°, 300°, and 330° can be analyzed using symmetry. With these settings, our study considered 8 sizes and 6 directions, leading to the construction of a total of 48 finite element models.

4.3. Numerical Calculations and Results Analysis

4.3.1. Material Parameters and Boundary Conditions

In the finite element models of rock mass analysis units, two sets of nodes and elements were defined: one set consisted of equal-thickness joint elements used for stress analysis of the fracture units, and the other set consisted of line elements used for seepage analysis of the fracture units. The mesh division for these units is depicted in Figure 8. The material constitutive relationship for the rock was assumed to be elastic, while the fracture unit followed the Goodman joint element constitutive model. Relevant material mechanical parameters are presented in

Table 2. Table for parameters of material property of the model.

Element Type	Deformation Modulus (GPa)	Poisson’s Ratio	Normal Stiffness (MPa/mm)	Shear Stiffness (MPa/mm)	Cohesion (MPa)	Friction Angle (°)	Density (kg/m3)
Rock	36	0.2	/	/	1.75	28.5	2.5 × 103
Joint	/	/	15	5.28	/	/	/

.

For the boundary conditions, as depicted in Figure 9, the model’s bottom was constrained as a fixed boundary, the left and right sides were subjected to horizontal link constraint boundaries, and normal loads ${\sigma }_y=1, 2, 3, 4 \mathrm{MPa}$ were applied to the top and bottom boundaries. Regarding seepage boundary conditions, the left boundary of the model had fixed water pressure $p_1=100 \mathrm{kPa}$, and the right boundary had fixed water pressure $p_2=0 \mathrm{kPa}$. The top and bottom were treated as impermeable boundaries.

4.3.2. Determination of Equivalent Deformation Parameters

Assuming that the rock is an impermeable medium and applying the loading and hydraulic boundary conditions to the rock mass analysis units, the seepage field within the fracture network is simulated. Taking III_2A rock mass under a load of ${\sigma }_y=2 \mathrm{MPa}$ as an example, the distribution of fracture water pressure in the rock mass analysis units on the 10th, 20th, and 50th days is shown in Figure 10a–c. From the figures, it can be observed that the fluid infiltration is very slow, and, with time, the pore water pressure gradually increases from left to right.

To study the deformation characteristics of rock mass analysis units of different sizes, a comparative analysis was conducted on the influence of different loads on the deformation distribution pattern. Using models with sizes of 6 m × 6 m, 10 m × 10 m, and 14 m × 14 m as examples, their deformation contours under load ${\sigma }_y=2 \mathrm{MPa}$ are shown in Figure 11. The observations include the following: ① in the 6 m × 6 m-sized model, the presence of fracture structural surfaces significantly influences the deformation distribution in the rock mass analysis units, resulting in notable discontinuities in displacement values on either side of the fracture structural surfaces. ② In the 10 m × 10 m-sized model, the influence of fracture structural surfaces on the deformation in the rock mass analysis units decreases, and the differences in displacement values on either side of the fracture structural surfaces also reduce. ③ In the 14 m × 14 m-sized model, the deformation in the rock mass analysis units exhibits a more uniform distribution, and the influence of fracture structural surfaces on the deformation is minimal.

By applying loading boundary conditions in various directions and sizes to rock mass analysis units, a simulation of the deformation field was carried out. The results, which are shown in

Table 3. Comparison of calculation results of equivalent deformation modulus of different rock mass analysis unit sizes (direction: 0°).

Size	4 m × 4 m	6 m × 6 m	8 m × 8 m	10 m × 10 m	12 m × 12 m	14 m × 14 m	16 m × 16 m	18 m × 18 m
Ex	3.6	4.4	5.2	5.0	5.4	5.7	5.6	5.5
Ey	2.2	2.7	3.2	4.0	4.4	4.5	4.8	4.8
Gxy	1.5	1.6	2.0	2.2	2.4	2.4	2.6	2.5

, include the variation curves of the x-direction-equivalent elastic modulus (E_x), the y-direction-equivalent elastic modulus (E_y), and the xy-plane-equivalent shear modulus (G_xy) with respect to rock mass analysis unit size, as depicted in Figure 12. The upper and lower limits of the effective equivalent parameters were established based on the calculation results for an 18 m × 18 m-sized rock mass analysis unit, using an acceptable error of 10% [38,39]. When the rock mass analysis unit size is taken as 10 m × 10 m or larger, the calculated values for E_x, E_y, and G_xy are all within the acceptable error range, and there are no significant fluctuations with changing rock mass analysis unit size, satisfying quasi-constant properties. In order to investigate the tensor-like properties, it is necessary to analyze the coefficients of variation in different directions (0°, 30°, 60°, 90°, 120°, 150°, 180°, 210°, 270°, 300°, and 330°) of the analysis units. The relationship between the coefficient of variation (Cv_x, Cv_y, and Cv_xy) and the rock mass analysis unit size is shown in Figure 13. When the rock mass analysis unit size is greater than 10 m × 10 m, all the coefficients of variation no longer exhibit significant changes, and the calculated values fall within the acceptable error range, satisfying quasi-tensor characteristics.

Following the aforementioned procedures and methods, calculations and analyses were conducted for several other types of rock masses, yielding the equivalent deformation parameters for various rock masses, as shown in

Table 4. Equivalent deformation parameters of various rock masses.

Rock Type	Rock Mass Rating	Weathering Degree	REV (m)	Ex (GPa)	Ey (GPa)	Gxy (GPa)
Potassium Feldspar Granite	IIA	Fresh, Slightly Weathered	10 × 10	15.7	15.3	7.1
	III1A	Mild Weathering (Lower Section)	10 × 10	7.2	5.9	2.8
	III2A	Mild Weathering (Upper Section)	8 × 8	5.3	4.6	2.1
	IV	Severe Weathering	8 × 8	1.67	1.43	0.7

. Due to limitations in experimental conditions, comparative tests on the deformation parameters of engineering rock samples were not conducted. However, the calculated parameters fall within the reasonable range of variations reported for granite deformation parameters in reference [40], validating the computational approach used in this study. It should be noted that physical and mechanical properties of rock masses vary across different regions due to differences in geological ages, weathering degrees, and degree of development of structural planes. Therefore, the calculated engineering rock mass deformation parameters in this study may exhibit some discrepancies when compared to those reported in other research.

4.4. Prediction Results of the Hybrid Model

Through numerical analysis, the coefficients for rainfall components in Equation (2) are determined by fitting the calculated results of the displacements of monitoring points under various rainfall conditions. The temperature- and time-dependent components are expressed using the statistical model expressions in Equations (3) and (4), and the model coefficients are obtained through stepwise regression analysis. By analyzing data on slope deformations from 2018 to 2022, the regression coefficients for the hybrid prediction model are obtained. Taking the horizontal displacements in the downslope direction at monitoring points M4xn1-1 and M4xn1-2 as an example, the coefficients are shown in

Table 5. Coefficients in the hybrid model of two monitoring points.

Coefficients	Monitoring Points		Coefficients	Monitoring Points
	M4xn1-1	M4xn1-2		M4xn1-1	M4xn1-2
a0	0.12373	−1.6799	b11	−0.17012	−0.29339
a1	0.06194	0.03898	b12	0.00689	−0.06974
a2	0.04831	0.05896	b21	−0.16921	0.05185
a3	0.02492	0.01244	b22	0.00728	0.15842
a3	0.00464	0.00651	c1	0.01342	0.02285
a3	0.00248	0.00793	c2	0.20623	0.41164

.

The fitting and predictive results are shown in Figure 14a,b. These figures illustrate that the predictive model proposed in this paper exhibits strong fitting and predictive capabilities, accurately reflecting the overall deformation patterns of the rock reservoir bank slope. However, the process of acquiring slope deformation monitoring data is inevitably affected by factors such as human error, external interference, and equipment failures. This results in irregular jumps in the monitoring data, particularly noticeable in the early stages of monitoring for measurement point M4xn1-2 (prior to 1 May 2019), leading to poor fitting performance during this period. Since there is no systematic pattern in the errors between measured and predicted data, no correlation with environmental variables, and no detected structural changes in the slope, these errors are likely caused by external random factors such as monitoring uncertainties and rock mass properties.

4.5. Comparison of Prediction Accuracies

To analyze the modeling effectiveness of this model, the complex correlation coefficient (R) and the root mean square error (RMSE) are introduced as evaluation metrics. The modeling effectiveness of the traditional statistical model M₁, the deterministic model M₂, and the hybrid model M₃ are compared, and the results are presented in

Table 6. Comparison results of different prediction models.

Evaluation Criteria	R			RMSE
Models	M1	M2	M3	M1	M2	M3
M4xn1-1	0.8941	0.8628	0.9358	0.1514	0.1961	0.1021
M4xn1-2	0.8077	0.8349	0.8627	0.4303	0.3439	0.2238

. It can be observed that the proposed hybrid model in this paper exhibits the highest R, closest to 1, and the smallest RMSE among all the models, indicating that this model provides the best prediction results.

Figure 15a,b presents a comparative diagram of the prediction results from the traditional statistical model, the deterministic model, and the hybrid model. The figures indicate that the proposed hybrid model achieves better prediction results than other comparative methods, demonstrating the effectiveness of combining numerical simulation with statistical regression. Figure 15c,d present the relative errors of prediction results from the three models. The figures reveal that, on certain dates (e.g., 8 February 2022, and 9 January 2023), the statistical or deterministic models perform better. However, for the majority of dates, the hybrid model demonstrates superior predictive capability. Additionally, due to the relatively modest measured average values (1.25 mm and 1.82 mm) during the forecasting period, coupled with substantial fluctuations (0.75 mm and 1.11 mm) and poor regularity, this results in notable relative errors.

Figure 16 displays boxplots of residuals from different models. It can be seen from the figure that the residuals at monitoring point M4xn1-1 are closer to zero, while those for monitoring point M4xn1-2 are more tightly clustered, with the median line of the residuals also approaching zero, indicating higher accuracy of this approach. It is noteworthy that all three prediction models exhibit certain biases. Referring to Figure 15, the authors suggest that this bias to the models mainly reflects the overall patterns in the training dataset, making them less effective in accounting for the occasional fluctuations in the validation dataset. For example, at the beginning of the validation set, without significant changes in environmental factors or structural conditions, M4xn1-1 abruptly increased from 1.39 mm to 1.59 mm, reaching the maximum value of the entire sequence; simultaneously, M4xn1-2 suddenly decreased from 2.42 mm to 1.62 mm.

5. Conclusions and Discussion

Considering the deformation characteristics of fractured rock masses under seepage conditions, this paper introduces the equivalent continuous medium theory. Through numerical simulations of discrete fracture networks considering the coupling of seepage and stress, the size effect and anisotropy of the deformation parameters were explored. Based on this analysis, a method for calculating the equivalent deformation parameters of fractured rock masses was proposed, consequently enhancing the predictive precision of the hybrid model. The principal conclusions can be summarized as follows:

(1) A method for iteratively solving the discrete fracture network model while considering seepage–stress coupling was proposed. The calculations of equivalent deformation and seepage parameters of the equivalent continuous medium under various load conditions were investigated, thus calculating the equivalent deformation parameters of the fractured rock mass;

(2) By harnessing the advantages of integrating statistical models with finite element analysis models, the proposed model can adeptly incorporate temperature- and time-dependent characteristics while also reflecting the mechanical mechanisms of rainfall factors on predictive quantities. Engineering case study results demonstrate that, compared to traditional statistical models and deterministic models, the hybrid model proposed in this paper exhibits superior predictive performance.

However, further elucidation is required regarding some assumptions and limitations of the proposed model. The deterministic component of this model only computes deformations within the elastic range of slope lines under rainfall effects, while attributing deformations caused by rainfall-induced groundwater level changes and other factors such as degradation and creep of slope rock masses to time-dependent effects. Whether the selected time-dependent components comprehensively reflect these influences remains to be studied. Additionally, slope monitoring often encounters sudden changes attributed to random factors. Currently, prediction intervals (PIs) are commonly employed to address this issue, which also remains a topic for future research.