Introduction

Over the last decades, soil nailing still attracts the attention of many engineers and researchers. It is a reinforcement technique that uses passive elements drilled and grouted sub-horizontally into the ground. This method is passive because it only resists when soil displacement occurs; it is also one of the techniques the most used in situ for soil retention around the world (Singh 2020). In addition, it is a practical technique used to stabilize the slope of vertical or nearly vertical cut, support excavations in the soil or in soft and weathered rock, the sidewalls of the access road for subways, the construction of tunnel gates, the construction of highways and railways, the abutments of bridges, etc.; finally, it is an innovative and very cost-effective technique.

The basic procedure for installing a soil nail includes excavation, drilling a hole, placing a nail in the hole, grouting the hole, and finally performing a facing operation. The working sequence is done from top to bottom.

In this paper, the reinforcement method of soil nailing walls was investigated, and the reinforcement parameters were studied. The effects of different reinforcement measures on the safety factor and cost were analyzed. Indeed, the determination of the optimal combination parameters with minimizing the objective function is an important problem of optimization; the goal is to achieve an optimal design that is safe, robust, and cost-effective. Recently, computing intelligence techniques such as artificial neural networks (ANN), genetic algorithms (GA), particle swarm optimization (PSO), ant colony (ACO), and design of experiments (DOE) such as Taguchi method, response surface methodology (RSM), and factorial design have been widely applied in many geotechnical optimization problems.

Sivakumar Babu and Basha (2008) introduced a target reliability approach for design optimization of concrete cantilever retaining walls; an optimization design method of composite soil nailing in loess excavation was introduced by Chang (2009); Ahmadi-Nedushan and Varaee (2009) proposed an optimization algorithm based on the particle swarm optimization (PSO) for optimum design of reinforced concrete Earth-retaining walls; Ghazavi and Bonab (2011) applied a methodology to arrive at the optimal design of concrete retaining wall using the ant colony optimization (ACO); Deb and Dhar (2011) proposed a finite difference-based simulation model and an evolutionary multiobjective optimization model to identify the optimal parameters for granular bed-stone column improved soft soil; Telis et al. (2013) presented a method for simultaneous optimization of the design characteristics of an Earth-retaining structure design using quality tools; Seo et al. (2014) estimated three design variables from the optimization design procedure proposed in soil nailing study considering constrained conditions; Zhang and Goh (2015 ) used a regression models for estimating ultimate and serviceability limit states of underground rock caverns; Yuan et al. (2019) proposed a statistical characterization of model which uncertainty has been exclusively focused on reliability-based design of soil nail wall; for reinforced concrete retaining wall design, Kalemci and Ikizler (2020) proposed an optimization algorithm called Rao-3; the total weight of the steel and concrete, which were used for constructing the retaining wall, were chosen as the objective function; Kalemci et al. (2020) used Grey wolf optimization algorithm of reinforced concrete cantilever retaining wall design.

In this study, response surface methodology (RSM) of design of experiment (DOE) was chosen as an optimization method; it allows mathematically rigorous and reliable research. A limited number of geotechnical studies have been reported in the literature on the use of RSM technique of optimization, from which we can cite the precursor work of Sivakumar Babu and Singh (2009) where, to better understand the deformation behavior of the soil nail walls, they carried out numerical experiments and regression models for the maximum lateral deformation and the prediction of the safety factor. Liu et al. (2013) proposed a response surface optimization design method for the foundation pit soil nailing, and Algin (2016) used RSM in the multiobjective optimization analysis to present the optimized design solution of a jet grouted raft. An interesting study by Güllü and Fedakar (2017) used RSM in the optimization of material stabilization in terms of their dosage rates. In their approach, Khoshnevisan et al. (2017) used response surface-based robust geotechnical design of supported excavation by means of a spreadsheet-based solution. Lafifi et al. (2019), on the other hand, used a combination of Taguchi’s design of experiment (DOE), RSM, and desirability function to optimize geotechnical parameters of a synthetic soil for pressuremeter test. In a recent research, Zhan et al. (2021) applied RSM to model unconfined compressive strength (UCS) of micaceous soils as a function of the dosage of various additives, an extensive study that has shown very interesting results. Sharma and Ramkrishnan (2020) presented a parametric optimization and multi-regression analysis for soil nailing using numerical approaches. Johari et al. (2020) used stochastic analysis using the random finite element method of system reliability analysis of soil nail wall.

The article is divided mainly into two parts; the first is devoted to the finite element model for representing the geometry and soil-structure aspects of a soil nailing system, while in part two is presented the optimization approach in which different parameters with their influence on the safety factor are discussed.

Most of the methods for soil nail wall stability analysis are based on the limit equilibrium method (LEM) for its simplicity and the reduced number of required parameters. A numerical study on the optimum layout of soil nails to stabilize the slopes was performed by Fan and Luo (2008); Tan et al. (2005) studied in detail the effects of 2D modeling of 3D soil nailing problem; Rawat and Gupta (2016a) used limit equilibrium and finite element methods in the analysis of a nailed soil slope; Deng et al. (2017) and Quansah et al. (2018) presented the optimum design with limit equilibrium analysis for the stability of soil nailed slope.

While the use of the finite difference or finite element methods has also attracted engineers in recent times, numerous researches have focused on FEM techniques for soil nail wall analysis; it is a powerful tool, which can be utilized for soil nail wall modeling. The main advantage of FEM is providing information about deformations of the soil nail system (Johari et al.2020).

Zhang et al. (1999) developed a three-dimensional (3D) finite element model (FEM) for the deformation analysis of nailed soil structures; Babu et al. (2007) analyzed a real vertical cut supported with retaining wall and soil nailing system, by 2D FEM program; Rawat and Gupta (2016b) analyzed using finite element software PLAXIS 3D the failure pattern and load-settlement plots for the various unreinforced and reinforced slopes. With the recent advances, the validation of soil nail models can be done using finite element software packages like Plaxis, Flac, GEO5, and Abaqus.

Material and methodology of analysis

Simulation model

In the present study, the finite element software Plaxis 2D was applied to perform the soil nailing simulation; it offers a number of advantages such as flexibility and operability. The geometry of the soil nailing system is presented in Fig. 1 with the parameters of interest.

Fig. 1
figure 1

Geometric model

The parameters included in the design process are the height of the soil nail wall H (m), the slope angle of the terrain β = 0°, the inclination of the wall α = 90°, the horizontal spacing of the nails Sh (m), the hole diameter D = 100 (mm), and nail diameter d = 32 (mm) in terms of (axial stiffness EA, bending stiffness EI).

In addition, the geomechanical data comprise the soil properties, such as the unit weight γ (kN/m3), the angle of the internal friction φ (°), the cohesion of the soil c (kPa), the dilatancy angle ψ (°), Young’s modulus E (kN/m2), and Poisson’s ratio ν (Fig. 2).

Fig. 2
figure 2

Finite element model for soil nailing

In the geometric model, fifteen (15) node triangular elements are used for generating a finite element mesh. Coarse mesh density was considered globally for the analysis but refined to a fine density mesh in areas surrounding each nail. Soil nail structures are modeled as plane strain problems in Plaxis 2D; using the elastic-plastic Mohr–Coulomb model, long-term behavior is simulated using drained analysis conditions.

To simulate a relatively larger area, the end boundary conditions were considered to be at roller condition on the sides (Ux = 0) and fixed at the bottom (Ux = Uy = 0).

In any soil nail wall modeling, staged excavation is needed and executed in a number of phases to simulate the construction sequences of the soil nail wall. The top-down construction sequence was simulated in the calculation stage using the staged construction technique available in Plaxis 2D.

Material properties

The most important input material parameters for plate elements are the flexural rigidity EI (bending stiffness) and the axial stiffness EA. The plate structural element is rectangular in shape with a width equal to 1 m in an out-of-plane direction.

The soil properties assumed in this analysis and the parameter of plate elements were used to model the nails and the shotcrete which are essentially based on Singh and Babu (2010).

The physical and mechanical properties of soil considered are given in Table 1.

Table 1 Soil model parameters (Singh and Babu 2010)

The nails and the facing of the soil nail wall were modeled using a plate element. Their stiffness was derived considering the grouted nails. Since the soil nails are circular in cross-section and placed at designed horizontal spacing, it is necessary to determine equivalent axial and bending stiffness for the correct simulation of circular soil nails as rectangular plate elements.

For the grouted nails, equivalent modulus of elasticity Eeq shall be determined accounting for the contribution of elastic stiffness of both grouts cover as well as reinforcement bar. From the fundamentals of the strength of materials, equivalent modulus of elasticity of grouted soil nail (Eeq) can be determined as (Babu and Singh 2009)

$$ {E}_{eq=\kern0.5em }{E}_n\left(\frac{A_n}{A}\right)+{E}_g\left(\frac{A_g}{A}\right) $$
(1)

where Eg is the modulus of elasticity of grout; En is the modulus of elasticity of nail; Eeq is the equivalent modulus of elasticity of grouted soil nail; \( A=\frac{\pi\ {D}^2}{4} \), where A is the cross-section area of grouted nail and D is the diameter of drill hole; \( {A}_n=\frac{\pi\ {d}^2}{4} \), where An is the cross-section area of nail and d is the diameter of the nail; Ag = AAn, where Ag is the cross-section of grout cover; axial stiffness EA [kN/m] = \( \frac{E_{eq}}{S_h}\ \left(\frac{\pi\ {D}^2}{4}\right) \); bending stiffness EI [kNm2/m] = \( \frac{E_{eq}}{S_h}\ \left(\frac{\pi\ {D}^2}{64}\right) \), where Sh is the horizontal spacing of soil nails.

In addition, equivalent plate thickness in meter is determined automatically by Plaxis 2D, \( {d}_{eq}=\sqrt{12\ \left(\frac{E_I}{E_A}\right)} \).

Nail and shotcrete properties taken as plate elements are given in Table 2, which have been calculated using the formulas given in Eq. (1).

Table 2 Material properties input for modeling

Safety factor for simulations was calculated using the strength reduction technique, i.e., a phi-reduction technique in Plaxis 2D, which in iterative steps reduces the soil shear strength parameters until the soil body fails (Plaxis Reference Manual 2002).

$$ \sum Msf=\frac{{\mathit{\tan}\varnothing}_{\mathrm{input}}}{{\mathit{\tan}\varnothing}_{\mathrm{reduced}}}=\frac{C_{\mathrm{input}}}{C_{\mathrm{reduced}}} $$
(2)

where ∅input is the input value of angle of internal friction (°), ∅reduced is the reduced value of angle of internal friction at failure (°), Cinput is the input value of cohesion (kPa), Creduced is the reduced value of cohesion at failure (kPa). Safety factors values are presented in Table 5.

Optimization model

Response surface methodology

Response surface methodology (RSM) was introduced by Box and Wilson in the early 1950s; it was a collection of statistical and mathematical techniques useful for the development and optimization of processes (Dwornicka and Pietraszek 2018). Various response surface methods are available; some of these are Box–Behnken design (BBD) and central composite design.

BBD takes the midpoints of the edges of the process space and the center point into consideration while constructing the design (Fig. 3); each variable was coded at the levels −1, 0, and 1, and it is composed of two parts: the central point and the middle points of the edges (Ranade and Thiagarajan 2017). This is classified as design points called space type. Therefore, the high value of the original variable is represented by (+1), and the low value is represented by (−1). The (0) value represents the average of these two values (Bagaber and Yusoff 2018). BBD is renowned as one of the most common and efficient designs used in the RSM. Using BBD facilitates the optimization of the effective parameters with a minimum number of experiments, as well as enabling analysis of the interaction between parameters (Polat and Sayan 2019).

Fig. 3
figure 3

Three-factor Box–Behnken design

Concepts and techniques of response surface methodology (RSM) have been applied successfully and extensively in different scientific and research areas (Zangeneh et al. 2002).

In this paper, Box–Behnken design (BBD) was applied because it is the most used and popularly RSM technique. It is more economical as compared to other designs due to the reduced number of experimental trials in the design.

The different steps involved in Box–Behnken design are given in Fig. 4.

Fig. 4
figure 4

Steps involved in Box–Behnken design

Parametric optimization

Several studies have been undertaken detailing the effect of parametric variations on the stability of nail slopes such as nail length, spacing between nails, nail inclination, and injection pressure. Fan and Luo (2008) investigated the effect of nail inclinations on the stability of soil nail slope by nonlinear FEM; Zhou et al. (2011) developed a three-dimensional (3D) finite element (FE) model to simulate the pullout behavior of a soil nail in a soil nail pullout box under different overburden and grouting pressures; Hong et al. (2017) conducted a total of eight soil nail pullout tests in a field to examine the influence of overburden pressure (OP) and grouting pressure (GP) on the PPR of soil nails. Jelusic and Zlender (2013) presented the optimization of ground nailing using a nonlinear programming approach to predict the safety factor using parameters such as slope angle, nail inclination, and length of nails. In addition, Seo et al. (2014) considered the optimization of ground nailing using three variable design parameters: bonded length of the nail, the number of nails, and the prestress; Sharma and Ramkrishnan (2020) undertook potential parametric optimization in soil nailing by considering soil nail interaction and back analysis of the pullout strength of nail using finite element analysis (Plaxis 2D).

Johari et al. (2020) used the model of Singh and Babu 2010 and represented the importance of the mentioned stability modes and their effects on system reliability and system probability of failure using the random finite element method.

In view of the above, the length of the nails, the vertical spacing between them, and the inclination of the nails were considered as the optimization parameters for the current study.

The process was carried out using Minitab 18 software (Minitab 2017) for selected factors to solve the optimization problem. The levels and ranges of factors used in this study are presented in Table 3 according to the guidelines (Clouterre 1991).

Table 3 Parameters and ranges for soil nailing analysis

One of the main objectives of this work was to determine the optimum values of variables for soil nailing wall. In this work, the desired goal in terms of stability and cost was defined as “minimize” to achieve the best optimum soil nailing parameters.

The results of the calculation of the safety factors were obtained based on the simulations carried out using finite element analysis with Plaxis 2D.

Regarding the calculation of the cost of the project, the total cost is considered as the sum of the cost of the various stages of construction (Table 4).

Table 4 Construction soil nailing steps

The total cost of the project is the total cost of construction work from top to bottom 1–6; in this case, the total cost is calculated without considering direct and indirect costs; the focus is on the cost optimization in terms of the installation of soil nailing alone. It should be noted that the cost estimation models adopted are only used for illustration purposes.

The cost of excavation (Cexc) per m3 is estimated at 400.00 DZD. For the 100-mm diameter drilling, the cost is 8000.00 DZD per meter of drilling (Cdrill). The diameter of the nail is 32 mm, and the cost of 1 m is 1500.00 DZD (Cnail). The cement grouting is estimated at 55,000.00 DZD per meter (Ccement). Cement grouting is necessary after drilling and installing nails to strengthen and protect them from rust. Shotcrete is estimated at 60,000.00 DZD per m3 (Cshot). Finally, accessories like the nail head are installed to provide protection to the tip of the frame.

$$ {COST}_{\mathrm{total}}={C}_{\mathrm{exc}}+{C}_{\mathrm{drill}}+{C}_{\mathrm{nail}}+{C}_{\mathrm{cement}}+{C}_{\mathrm{shot}} $$

The total cost shown below (Table 5) is expressed in Algerian Dinar (DZD).

Table 5 Experimental results according to Box–Behnken design

Results and discussion

In this study, BBD was carried out with a total of 15 experiments. Based on the Tagushi method, we could have used 27 experiments (L27) but thanks to Box–Benhken design the number was reduced to 15 experiments with very satisfactory results.

The results are analyzed so that the conditions could be optimized to give the best outcome for all response parameters. Results of calculation of safety factors and total cost are summarized in Table 5.

Development of a regression model equation

The experimental data of independent variables and responses were fitted in a 2nd-order polynomial equation which is given in Eq. (3).

$$ Y={B}_0+{\sum}_{i=1}^K{B}_i{x}_i+{\sum}_{i=1}^K{B}_{ij}\ {x}_i^2+{\sum}_{i=1}^{K-1}{\sum}_{j=2}^K{B}_{ij}{x}_i{x}_j $$
(3)

where Y is the predicted response (safety factor) or (cost), B0, Bi, Bii, and Bij are assigned to constant, linear, quadratic, and interactive values of regression coefficients, respectively, Xi and Xj are input variables, and i and j are the index numbers.

The developed regression equation represents the quantitative effect of input factors and their interactions.

Based on the experimental design given in Table 5, the RSM proposed Eqs. (4) and (5) are mathematical models to describe the effect of investigated factors on safety factor and cost, respectively.

$$ Fs=1.13+0.072\ L-0.0062\ i+0.139\ {V}_s-0.00262\ {L}^2-0.00052\ {i}^2-0.101\ {V}_S^2\kern0.5em +0.00207\ L\ast i+0.0044\ L\ast {V}_s+0.0032\ i\ast {V}_s $$
(4)
$$ Cost=-934753-373202\ L+454581\ i+3645800\ {S}_V+43709\ {L}^2-6994\ {i}^2-399351\ {S}_V^2-75000\ L\times {S}_V-139871\ i\times {S}_V $$
(5)

It was found that the safety factor Eq. (4) had an R2 correlation coefficient value of 0.90, and the cost model Eq. (5) had an R2 value of 0.94. It can be verified from the R2 values that the predictions of the model equation corresponded to the experimental values.

According to Table 6, a good agreement between the experimental and predicted values confirms the validity of the model.

Table 6 Experimental and predicted results according to Box–Behnken design

Analysis of variance

Analysis of variance (ANOVA) is a statistical technique carried out to find the most significant process parameters which will affect output parameters. The ANOVA computes the quantities such as degrees of freedom (DF), multiple determination coefficients (R2) tests, F ratio (F), and p-value (p), as illustrated in Tables 7 and 8, which were used to determine the model’s significance.

Table 7 ANOVA for response surface quadratic model (safety factor)
Table 8 ANOVA for response surface quadratic model (cost)

If the p-value is less than 0.05, the process parameter is said to be significant, and if the p-value is bigger than 0.05, the process parameter is said to be insignificant.

The results of ANOVA presented in Tables 7 and 8 reveal that the most effective parameter is the length, as depicted by its low p-value (p < 0.05) and high F value, then followed by vertical spacing with relatively high p-value and considerably low F value, and finally nail inclination factor.

The suitability of the quadratic response surface model was justified by ANOVA. Only the length of nails parameter has a significant influence on the safety factor (Fig. 5a). In Fig. 5b, we can see that length and vertical spacing have a main effect on cost.

Fig. 5
figure 5

a Main effect plot for safety factor. b Main effect plot for cost

To determine the magnitude, direction, and importance of the effects, we can use the normal probability diagram of the effects. On a plot in Fig. 6a, the main effect for factor A (length) on the safety factor is statistically significant at the 0.05 level. This point has a distinct color and shape (red square) from the points for the insignificant effects (blue circle). Whereas on plot Fig. 6b, the main effects for factors A (length) and C (vertical spacing) on cost are statistically significant. The normal plot of the standardized effects in Fig. 6a and Fig. 6b proved that the most significant factor on safety factor is length, and the most significant factors on cost are the length and vertical spacing.

Fig. 6
figure 6

a Normal plot for Fs. b Normal plot for cost

The Pareto chart was developed to compare the relative magnitude of the effects of various factors on the response, their significance and interactions. The effects are plotted in decreasing order of the absolute value of standardized effects and draws a reference line on the chart (Dwornicka and Pietraszek 2018).

As shown in Fig. 7a, the input effect Pareto bar A (length of nails) is to the right of the vertical red line; therefore, this bar is statistically significant on safety factor at a 5% significance level with the current model terms. Although, the other factors C, CC, AB, BB, B, AA, BC, and AC seem insignificant. On the other hand, in Fig.7b, A and C bars are statistically significant on the cost response variable.

Fig. 7
figure 7

a Pareto chart (response variable Fs). b Pareto chart (response variable cost)

Optimization using response contour lines

The relationship between factors and responses is well understood when using contour lines plots. The contour plot provides a 2D view of the surface where points that have the same response are connected to produce contour lines of constant responses. In Fig. 8a–c, the contour plots show the response as a function of the other variables: length, inclination, and vertical spacing. The plots are very explicit regarding the influence of the parameters on the safety factor. We can clearly see the increase in the safety factor with the increase in length and lessened inclination. On the other hand, the safety factor decreases with the increase in vertical spacing (Fig. 8c). Figure 9a–c plotted the second response variable (cost) in combination with the rest of the parameters. The cost decreases with a decreasing length while vertical spacing is kept constant and inclination has little effect (Fig. 9a). Vertical spacing (Sv) affects the cost; the smaller Sv, the higher the cost as more nails will be needed (Fig. 9c). We can see that inclination has a little effect on the cost also.

Fig. 8
figure 8

a Contour lines (Fs) vs inclination and length. b Contour lines (Fs) vs length and vertical spacing. c Contour lines (Fs) vs inclination and vertical spacing

Fig. 9
figure 9

a Contour lines (cost) vs inclination and length. b Contour lines (cost) vs length and vertical spacing. c Contour lines (cost) vs inclination and vertical spacing

Using the response optimizer to identify the combination of input variable settings that optimize two responses, Fig. 10 shows that the values of response variables can be predicted and obtained from the developed model. An optimal combination was found with a length of 8 m, an inclination of 25°, and a vertical spacing of 2 m to allow a safety factor of 1.5721 and an estimated cost of 3.371 × 106 DZD/linear meter with good desirability of 0.9939.

Fig. 10
figure 10

Optimization

Once the optimal combination of the design parameters has been obtained, the next step is to predict and verify the required improvement using the optimal combination of design parameters (Table 9).

Table 9 Results of the confirmation experiment

Good agreement between the predicted and experiment safety factor and cost was observed.

Conclusion

The present study undertakes parametric optimization in soil nailing by considering soil nail length, vertical spacing, and nail inclination. The mechanical model is based on finite element analysis, and the optimization was carried out using the RSM method.

The study demonstrated conclusively that the use of the Box–Behnken design and RSM for the optimization of soil nailing parameters is a very useful statistic method that reduces the number of experiments.

From ANOVA data, it was shown that the change in the nail length has a more significant impact on the stability of the soil nailing and cost, followed by the vertical spacing and inclination. The safety factor increases with an increase in length due to an increase of axial nail force, shearing force, and bending moments to resist the loading and deformation. The cost is influenced by the vertical spacing of nails as well as their length.

The mutual interactions between the independent variables are described with quadratic equations which makes it possible to predict the response under the applied conditions. However, one should note that the response surface models obtained are valid only in the selected ranges of the parameters.