Determining Groundwater Table in Slopes With Horizontal Drains and Design Framework for Length and Spacing

Horizontal drains (HDs) were initially introduced for groundwater management in agricultural land (Stuyt et al., 2005). Today, they are widely used as an effective technique for drainage management in a variety of geotechnical applications, including slope and embankment stability improvement (Nonveiller, 1981; Lau and Kenney, 1984; Martin et al., 1995; Cai et al., 1998; Forrester, 2001; Rahardjo et al., 2003, 2011; Tsao et al., 2005; Ghiassian and Ghareh, 2008; Shrestha et al., 2008; and Hinds et al., 2021), drainage management behind soil nailing and retaining walls (Chen and Chen, 2013; Geotechnical Engineering Office, 2017; Lakruwan and Kulathilaka, 2020, 2021), enhancement of the liquefaction resistance of soil (Fasano et al., 2019, 2021), and acceleration of consolidation (Feng et al., 2020). This study focused on evaluating the effectiveness of HDs in enhancing slope stability.

The effectiveness of HDs in enhancing slope stability is primarily based on several design parameters such as location, spacing, length, and drain type (Rahardjo et al., 2011; Lakruwan et al., 2021; and Abed et al., 2024). Rahardjo et al. (2003) showed that installation of HDs at the bottom of a slope was more effective than installation of HDs at higher elevations, thus confirming the results of Lau and Kenney (1984) and Martin et al. (1995). Additionally, some studies have tried to identify a maximum limit for the length of HDs (Royster, 1980; Lau and Kenney, 1984; Cai et al., 1998; and Cook et al., 2008), and several attempts have been made to develop design charts (Kenny et al., 1977; Stanić, 1984), but they are generally only valid for the specific slope conditions (Resnick and Znidarčić, 1990). Additionally, the general guidelines provided by Cook et al. (2012) for HD spacing based on the soil type vary significantly depending on the soil and recharge conditions.

Previous studies have proven that the design of HDs is highly site-specific and that providing generalized guides is not feasible. Hence, instead of a generalized guide, a novel method that is applicable to a wide range of slope conditions is proposed in this study to uniquely analyze the effects of the length and spacing of HDs on the pore water force reduction on a slip surface and, subsequently, the stability enhancement of a slope.

Although a few studies (Greg et al., 2013; Conte and Troncone, 2018) have used the Hooghoudt (1940) method for slope stability applications, its applicability to slopes is limited because the assumption of semi-infinitely long and parallel HDs is invalid for most slope stabilization applications. Although the Hooghoudt (1940) method assumes a uniform variation in the GWT profile along the HDs, the limited lengths result in non-uniform GWT variations along the HDs on slopes, leading to a complex distribution of the GWT profile in three-dimensional (3D) space.

Crenshaw and Santi (2004) were the first researchers to address this issue and proposed an empirical method for determining the 3D variation in the GWT along short HDs. However, the requirement of a GWT height of approximately 20 ft (6 m) behind the drain field and the inability to predict the GWT profile beyond the HD field are drawbacks of Crenshaw and Santi’s (2004) method.

The accumulation of pore water pressure (PWP) on the slip surface can be identified as the major triggering factor for rain-induced slope failures. Several studies have evaluated slope stability based on the pore water force on the slip surface (U) (Tsao et al., 2005; Conte and Troncone, 2012, 2018).

Currently available analytical and empirical methods do not adequately address the effect of 3D variation in the GWT in short HDs on slopes. Additionally, two-dimensional (2D) flow simulation models neglect the 3D effect of the GWT profile, potentially resulting in a inadequate design (Lakruwan et al., 2022; Abed et al., 2024). Consequently, to accurately capture the comprehensive effect of the 3D variation of the GWT, complex 3D flow simulation models are necessary. However, setting up an accurate 3D model to represent an actual slope is challenging, and the time and computational power required to iterate such complex models to obtain the optimal design parameters for HDs significantly increase the complexity of the process.

In contrast, this study proposes a novel semi-empirical method for predicting the steady-state GWT profile of HDs with limited lengths based on a series of numerical simulations. This method was validated through a series of numerical cases and field studies. The proposed method has fewer limitations, requires minimal computational power, can be conveniently implemented, and not only can predict the GWT profile in the HD field, but also the slope behind the HD.

Furthermore, we extend the proposed GWT prediction method to obtain the U and reduction of U (ΔU) by HDs. We then describe the procedure for developing a design chart to select the optimal length and spacing of the HDs for a given slope. To demonstrate the effectiveness of our proposed method and design procedure, we performed two case studies, and we provide several guidelines for selecting the optimal length and spacing of HDs.

The proposed GWT prediction method and design procedure allow designers to efficiently perform HD design while considering the 3D nature of the problem.

The objectives of this study are summarized as follows:

1. To develop a novel semi-empirical method to predict the GWT profile of a slope with HDs based on the observation of numerical simulations and to validate the proposed method through a comprehensive series of numerical case studies and a rigorous field study, ensuring the reliability and applicability of the proposed method;

2. To establish an average GWT profile that adequately captures 3D variations, which may be utilized in 2D numerical simulations, providing a more accurate and efficient approach for analyzing slope stability;

3. To quantify the U and ΔU values in the presence of a slip surface by integrating across the average GWT, and to thoroughly analyze the effect of the HD length and spacing on ΔU, highlighting the significance of selecting suitable values to achieve effective and economical design; and

4. To produce a design chart for specific slope conditions, enabling designers to confidently select optimal HD lengths and spacing based on the design requirements, and to provide insights and recommendations for the selection of HD lengths and spacing to achieve an optimal design.

The structure of this article (Figure 1) is summarized as follows:

1. Proposed Method: The novel semi-empirical method developed to determine or predict the GWT on a slope with HDs based on a series of 3D finite element (FE) analyses conducted on slope I is presented and shown in Figure 2. The average GWT obtained from the results is subsequently used to calculate the U and ΔU.

2. Validation: The proposed method for predicting GWT, U, and ΔU is comprehensively validated through a numerical case study on slopes I and II in Figure 2 and a field study.

3. Results and Discussion: The effect of the HD length and spacing on the ΔU (and factor of safety [FS] enhancement, ΔFS) value is rigorously evaluated using slopes I and II in Figure 2.

Series of FE steady-state flow analyses were performed on a 30-m-wide (in the Y direction) section of slope I using PLAXIS 3D (CONNECT Edition V20) (PLAXIS, 2020). Ten-node tetrahedral elements were used for FE discretization. The 3D FE model of slope I is shown in Figure 3. Table 1 presents the details of the numerical model analysis, and the parameters are shown in Figures 3 and 4.

Establishing the upstream hydraulic boundary in slopes is challenging due to its complex behavior influenced by hydrological conditions outside the slope domain. Following previous studies (Das et al., 2022; Lernia et al., 2022), a constant-head Dirichlet boundary with a total head of H_i (D + h_i, where h_i is the height of the GWT at the boundary relative to the HD level) was used. Although this upstream head exhibits seasonal variation (Lernia et al., 2022), using the high head value of the rainy season is a conservative approach in steady-state analysis. The constraints for the location of the upstream boundary are outlined in the section on “Constraints of the Proposed Method.” A Neumann boundary condition was assumed at the top surface to simulate the infiltration during rain. Although infiltration has temporal variability, we used the average value of the peak rainfall as a constant recharge (R) to provide a conservative estimate in steady-state conditions. This recharge represents the infiltration after surface runoff. Additionally, a 0.1 m cutoff value was used to prevent excessive ponding on the surface. Line drains with zero water pressure were employed as the HDs, neglecting the effect of entrance head loss and pressure head inside fully flowing pipes.

We observed three distinct regions of groundwater flow on slopes with HDs, which were characterized by flow direction (Figure 5). In the region behind the HDs, the X-directional component of the groundwater flow (q_x) dominates the Y-directional component (q_y), resulting in a 2D flow in the XZ plane. Conversely, when HDs were present, q_y became dominant over q_x, leading to 2D flow in the YZ plane. The transition from the XZ plane to the YZ plane occurred through a narrow region at the end of the HDs. Because clear distinctions between the three regions could not be established based solely on groundwater flow, we examined the GWT profile at the mid-section between the HDs (on the XZ plane).

The gradual increase in the GWT in the HD field was idealized into two linear profiles with different gradients to simplify the proposed method. Zone 1 corresponds to the region with flow in the YZ plane and exhibits a gradual increase in GWT with a small gradient. Zone 2 was located near the drain end, where a sudden increase in the GWT occurred, corresponding to the transition region. The length of this region is approximately equal to the drain spacing (S) and has a steeper gradient than that of zone 1. Finally, zone 3 represents the curved GWT profile behind the HD field, where the flow is in the XZ plane.

In summary, we defined three distinct zones in the GWT profile based on our observations of flow direction and changes in the GWT gradient along the profile. These zones are labeled as zones 1, 2, and 3, as illustrated in Figure 4.

The GWT profile in the mid-section between HDs can be predicted as follows.

• Zone 1:

where Z_w is the GWT profile (m), and D denotes the depth to the impermeable layer (m).

• Zone 2:

and

• Zone 3: Zone 3 may be characterized by two different methods, in local coordinates and global coordinates.

HDs are typically inclined to facilitate gravity drainage, resulting in a higher GWT profile than horizontally installed HDs. Therefore, we introduced a correction to the proposed method to account for the inclinations of the HDs.

In this GWT prediction method, zones 1 and 2 provide empirical solutions, whereas zone 3 represents a theoretical GWT profile with an empirical boundary (h_e).

The accuracy of the proposed method for predicting the GWT at the mid-section between the HDs was validated through a comprehensive series of numerical case studies using the 3D FE model of slope II and a field study.

The GWT profile variation between the HDs on the YZ plane for slope I having HDs with 10 m spacing and 30 m length is shown in Figure 7. The GWT profile between the drains conformed to the half-elliptical pattern proposed by Hooghoudt (1940) along the length of the drains. However, a uniform GWT profile was observed behind the HD field. Hence, the average value of the GWT profile along the length of the HDs (zones 1 and 2), denoted as (Z_w(avg)), as shown in Figure 8, can be used to characterize the GWT profile for further calculations of U, ΔU, and factor of safety (FS) enhancement (ΔFS), as well as incorporation into 2D numerical simulations (see Supplemental Material Appendix A).

The constraints of the proposed GWT prediction method and concomitant design chart are as follows:

• The slope is assumed to consist of isotropic and homogeneous soil.

• A horizontal impermeable layer is present at the bottom of the slope.

• A known groundwater level exists sufficiently away from the slope that can be used as the upstream constant head boundary. It is recommended to locate this boundary at least 15 m away from the HD field, which is the minimum distance considered in this study.

• An equivalent constant recharge rate must be defined for infiltration.

• There is a uniform distribution of water level in the Y direction. In the case of a non-uniform distribution, particularly h_i, it is recommended to divide the slope into sections with uniform water level and apply the method to each section separately.

The proposed method involves multiple rigorous iterative calculation steps to determine the optimal length and spacing of HDs. To overcome the time-consuming nature of the design process, a Python-based (Python 3.6.2; Rossum and Drake, 2009) computer program was developed to automate the calculations. The step-by-step design procedure for HDs, based on the framework proposed in this study, is outlined in Figure B1 (Supplemental Material Appendix B).

This section presents the verification of the proposed GWT prediction method using a series of numerical case studies and a field study. Also, the accuracy of the calculation method for U and ΔU is discussed.

Figure 10A shows a comparison of the predicted and 3D FE model GWT profiles for horizontally installed HDs with different lengths and spacings. The small MAE values indicate good agreement between the GWT profiles. Generally, zones 1 and 2 showed less agreement between the GWT profiles than zone 3, which can be attributed to the empirical linear idealization of the GWT in zones 1 and 2. In contrast, the theoretical prediction of the GWT profile with an empirical boundary in zone 3 resulted in lower MAE values. Furthermore, the agreement between the GWT profiles deteriorated with an increase in the spacing and length of the HDs. When the HD length increased, the distance to the upstream boundary of the slope (distance to h_i) decreased, leading to poor predictions. Therefore, it is crucial to establish the h_i at a sufficient distance away from the expected length of the HDs to ensure accurate results.

Figure 10B compares the predicted and 3D FE model GWT profiles for horizontally installed HDs at different rainfall intensities (n). Higher rainfall intensities resulted in a lower agreement between the GWT profiles, and a slight overestimation in zone 2 was also observed. The lower conformity indicated by the lower R² value (Figure 11) for the prediction of h_cr can be identified as the reason for this behavior. However, it did not significantly affect the accuracy of the GWT prediction. Figure 10c and d confirms the accuracy of the proposed method in predicting the GWT under different initial GWT conditions and HD angles, respectively.

Figure 11 compares the h_c and h_e values predicted by the proposed method with the values obtained from the 3D FE models. The low MAE and high R² values suggest that the proposed method provides a good prediction of h_c and h_c.

The Ginigathhena landslide, located at bridge no. 48/2 at the Avissawella-Hatton-Nuwaraeliya Road in Sri Lanka, was selected for field validation of the proposed method. The details of the landslides were documented by Lakruwan (2019) and Lakruwan and Kulathilaka (2020, 2021). The landslide mass consisted of thick colluvium and a residual soil layer underlain by highly crystalline metamorphic rock. The initial slope failure occurred on June 22, 2014, followed by excavation at the toe for road widening and propagation to the upper slope during the rainy season. The landslide was mitigated by improvements in drainage and soil nailing.

The average homogeneous and isotropic hydraulic conductivity of the soil overburden was found to be 8 × 10⁻⁶ m/s, and the bedrock was considered to be impermeable (Lakruwan, 2019). The water levels of the two observation wells, W1 and W2, located 1.7 m and 42.38 m away from HD field, respectively, were monitored after the 20-m-long HDs with 5 m spacing were installed (Figure 12). Because the horizontal impermeable bottom layer only existed in that section, only the top slope section demarcated by the dotted line rectangle in Figure 12 was considered for the validation of the proposed method. Two cases were selected for validation: Case 1, no rainfall from February 27 to March 21, 2018, and Case 2, 450 mm cumulative rainfall from May 18 to 23, 2018 (n = 0.11). Since the water level of W2 was rather constant with time, it was selected as the H_i (B = 62.38 m and H_i = 15.00 m [Case 1]/15.58 m [Case 2]). The observed and predicted water levels of W1 in Case 1 (zero elevation at the top of the impermeable layer) were 2.67 and 2.77 m, respectively. The observed and predicted water levels of W1 in Case 2 were 3.61 and 4.06 m, respectively. Therefore, despite a slight overestimation, the proposed method can accurately predict the water level.

Slope I was used for the validation of the proposed method for obtaining U and ΔU. The average U (kN/m) from the 3D FE model was obtained by integrating the PWP values of the 51 stress points located on the FE mesh of the slip surface (Figure B2 in Supplemental Material Appendix B). Three point-rows were located, one on a section on the HD and one each on sections at the mid- and quarter-points between the HDs. Each row consisted of 19 points at 2.5 m spacings in the X direction.

The good agreement between the U values obtained from the 3D FE models and the proposed method (Eq. 24 and Eq. 25) for slope I (h_i = 6 m and n = 0.1), as shown in Figure 13a, confirms its validity. Also, the plots of ΔU in Figure 13b confirm the accuracy of the proposed method in calculating ΔU.

This comprehensive validation confirmed the accuracy of the proposed GWT prediction method for slopes with HDs under a wide range of conditions subjected to the aforementioned limitations. Moreover, it validates the prediction of U and ΔU in the presence of predefined slip surface. Therefore, the proposed method can be used as an efficient tool for designing HDs while overcoming the difficulties of complex 3D numerical simulations.

This section presents examples of applications of the proposed method for the design of HDs and provides several conclusions on selecting appropriate lengths and spacing for HDs.

A 3D FE stress-deformation analysis was performed on slope I using the strength reduction method to obtain the PWP on the slip surface and FS. Figure 14 shows the PWP distribution on the slip surface obtained from the 3D FE model of slope I with a drain spacing of 10 m. The PWP decreased with HD length, and no further reduction occurred beyond normalized drain length (L/Ls, where Ls is the length to the slip surface) of 2.0. The effective stress increased, as shown in Figure 15, while complying with the decrease in PWP, which eventually resulted in similar behavior in the FS enhancement.

Figures 13 provides further insights into the effects of the HD length. In Figure 13, both U and ΔU achieve uniform values for long HDs, and a further increase in the drain length becomes insignificant, which confirms the results of Cai et al. (1998), Cook et al. (2008), Lau and Kenney (1984), and Royster (1980). This phenomenon is explained in Figure B3 (Supplemental Material Appendix B). In short HDs, the failure wedge extends into all three zones of the GWT, and U becomes highly sensitive to the length of the HDs because the high GWT in zone 3 gradually moves away from the failure wedge. With an increase in HD length, zones 3 and 2 gradually move away from the failure wedge, eventually limiting the failure wedge to zone 1. In zone 1, U shows negligible sensitivity to drain length increments because the gradient of the GWT is very small.

Figure 16 illustrates the effect of n, ranging from 0.01 to 0.2, on U for different drain lengths and spacing. The U values at n = 0.01 and n = 0.2 are denoted as U_0.01 and U_0.2, respectively. High GWTs in widely spaced HDs are highly sensitive to surface recharge, making n significantly affect U at large S values. This effect reduces as the GWT lowers in closely spaced HDs. Furthermore, the impact of n on U is significant in short HDs and decreases with increasing drain length, due to the shift in the slip surface location from zone 3 to zone 1, as discussed in the previous section.

A design chart can be used to efficiently select the length and spacing of HDs to achieve the desired ΔU or ΔFS values for a particular slope failure. The procedure for calculating ΔFS is outlined in Supplemental Material Appendix C. Given the heterogeneity of the failures, the design charts by Kenny et al. (1977) and Stanić (1984) were developed for unique slope conditions. Therefore, instead of producing generalized charts, we propose creating unique design charts for a given slope to select the drain length and spacing. The limitations on the chart used are outlined in the section “Constraints of the Proposed Method.”

Figure 17 shows the design charts, and variation in ΔU and ΔFS with drain spacing for different L/Ls, produced for slopes I and II. The non-compatibility of the design charts for slopes I and II shows that they are unique for a given slope, and a generalized design chart cannot be adopted.

A reduction in the sensitivity of the drain length for long HDs is also reflected in the design charts, where the design charts indicate that the effect of drain spacing is more significant in short HDs, whereas it is not significant in long HDs. In long HDs (i.e., when the failure wedge is located in zone 1), increasing the spacing causes only a slight increase in the GWT profile inside the failure wedge. Therefore, only a small increment in U is observed when increasing the spacing in long HDs, which is also shown by Abed et al. (2024).

A target value of ΔFS or ΔU can be achieved with a smaller HD length per unit slope width when large drain spacings are used (Figure 18). For example, to achieve ΔU of 9.9 kN/m² (ΔFS ≈ 0.6) in slope I, 5.3 and 1.7 m pipe lengths per meter width of slope are required when the spacing is 1.0 and 5.0 m, respectively. Therefore, the use of a larger drain spacing is more economical. However, two main issues arise with a wide drain spacing. One issue is the risk of non-functioning HDs for various reasons, such as clogging of the filter screen, mineral precipitation, and root growth inside the pipe (Martin et al., 1995; Lee, 2013). In the case of a non-functioning drain pipe, the drain spacing will be doubled, and ΔFS will be reduced to 0.54 and 0.35 for the initial drain spacing of 1.0 and 5.0 m, respectively. Therefore, a small spacing is more conservative in such situations. This issue can be eliminated using long HDs because the effect of spacing is not significant. The next issue associated with widely spaced HDs is the requirement for a pipe with a higher discharge capacity. Therefore, a trade-off exists between the pipe flow capacity and total pipe length, and the flow capacity of the drain pipe should also be considered in the selection of drain spacing.

The following conclusions have been drawn from this study.

1. A novel semi-empirical method that can be used to obtain the GWT profile of a slope with HDs subjected to constant recharge from the surface is proposed in this study, along with a design framework for HD length and spacing. The accuracy of the proposed method was comprehensively validated through a series of numerical case studies for a wide range of conceivable conditions and monitoring data of an actual slope.

2. The proposed method was extended to determine the PWP force on a predefined slip surface via linear integration between Z_w(avg) and Z_s, which was subsequently verified through a numerical case study. The Z_w(avg) can be used in 2D numerical models to replicate the 3D effects of the GWT.

3. A Python-based computer program was developed to automate multiple rigorous iterative calculation steps for determining the optimal length and spacing of the HDs. A unique design chart was produced for a specific slope using the developed computer program.

4. Design charts of slopes I and II were used to evaluate the effects of length and spacing of the HDs on the PWP reduction and FS increment of a slope. The HDs must only be extended to a length equal to the drain spacing beyond the slip surface. The length and spacing of the HDs required to achieve the target FS increase can be obtained from the design chart produced for a particular slope condition. Although lengthening the HDs is more effective than reducing the drain spacing, the risk of drain blockage and discharge capacity of the drain pipe must be considered when selecting the drain spacing.

Supplemental Material associated with this article can be found online at https://www.aegweb.org/e-eg-supplements.

This research was supported by the Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for Scientific Research (KAKENHI) grant (no. JP24K00974).