Numerical study for optimal design of soil nailed embankment slopes

This paper presents a set of new required programmes written in format data using the built-in programming FISH language available in three-dimensional Fast Lagrangian Analysis of Continua software (FLAC3D). These script data files were developed, to overcome difficulties noted during the nailed slope stability prediction. They established to analyse stability of general embankment slope cases under various soil nailing parameters design. To deal with 3D finite difference analysis results in terms of the factor of safety (FOS) and critical slip surface, both 2D finite difference and limit equilibrium (LEM) methods are also conducted. A comparison of the 3D and 2D results indicates that 2D analyses underestimate seriously stability predictions, notably in terms of the FOS. Therefore, it was concluded that, 3D analysis should be used as an alternative to the 2D analyses, even though it requires more time. Some useful conclusions are obtained from this work, which provides a good practice guidance for real life scenarios.

geotechnical engineering practice [9][10][11]. In all of the above studies, the design of these reinforced structures are usually limited to simplified cases with few rows and columns of soil nails. However, the behaviour assessment of soil nailed slopes with multiple rows and columns, that represents the general case, has received a very little attention. Moreover, the difficulty to create required FLAC 3D programmes to suit specific needs for general cases hinders seriously the application of numerical modelling in soil nailing technique. Therefore, the lack of these suitable programmes of soil nailing process requires further studies.
The main objective of this paper is therefore to illustrate the soil nailing design effects on the stability of an embankment slope resting partially on soft clay layer, having a similar geometry that studied by Cala et al. [8]. In fact, two different inclusion types that contains the embankment slope, namely one is working for the shear strength improvement (soil nails) but the second acts against it (soft clay layer) were considered in this study. On the other hand, both SLOPE/W and FLAC 2D were used to deal FLAC 3D results. The comparison is mainly established in terms of the FOS and critical slip surfaces. In order to manipulate different soil nailing design parameters such as, nail inclination angle, nail head size and soil nail layout, several new specific programmes using the programming FISH language available in FLAC 3D software, were developed by authors in this paper. The layout of the corresponding geometric model used in the present study is shown in Fig. 1.

Geometric model and soil properties
The geometric model of the soil nailed embankment slope considered in this study is 10 m height, with 45 • as a slope angle. First, nails of 12 m length are installed at 1.5 and 2 m spacing in horizontal and vertical directions respectively. The diameter of the used steel bar is 40 mm while, the grout hole diameter is 80 mm . The embankment slope rests partially on a thin soft clay layer of 1 m thick and 14 m depth (in x-direction). The water level was assumed to be at the ground surface. Due to the problem symmetry, only one-half of the model is simulated. Although, the embankment fill consists of gravely clay, a marly clay layer is considered as the main embankment foundation soil as shown in Fig. 1. Since these structures generate positive excess pore pressures within the foundation soils, the most critical condition occurs during the construction stage. Therefore, the model was run under fully undrained conditions, during which, soils are modeled with the elastic perfectly plastic Mohr-Coulomb (M-C) constitutive model, using drained parameters. Soil layers and their properties are shown in Fig. 1 with values provided in Table 1.

Limit equilibrium modelling
The computer software package SLOPE/W developed by GEOSLOPE has been used in this study. SLOPE/W packages enables construction of the soil embankment slope by defining its regions. Once the geometry of regions have been outlined, soil properties are assigned. Phreatic level is then added to numerical model via the PWP tab, available from the Analysis Settings in the KeyIn menu [12]. Subsequently, the numerical model is analysed under different combined length and inclination angle nails. However, the nail lengths used in this study are 9, 12, 15, and 18 m . In turn, every nail length is still combined with seven different nail inclination angles from horizontal of, 0, 5, 10, 20, 30, 40 and 60 • respectively. In all, 28 soil nailed embankment slope models are established and analysed in order to attempt an optimal nail inclination angle. Furthermore, in order to perform the nail layout analysis, eight other soil nailed embankment slope models are established. In fact, four nail inclination angles of 0, 5, 20 and 30 • are used respectively. In the first case, the shorter nail is installed at the upper part of the embankment slope, while the longer at its bottom. In which, the nail lengths from the top to bottom are 6, 9, 12, 15 and 18 m . In the second case, the design is reversed, so the shorter nail is at the bottom but the longer is at the top. Moreover, yes or no anchorage options are available as tools in SLOPE/W to allow simulation of fixed or free nail head modes. The nail forces are modeled as distributed forces over the nail length and the overall global factor of safety FOS is included in the analysis using F of S dependent option. It means that the applied nail load is actually the mobilised reinforcement load and is included in the FOS computations in the same manner as the soil mobilised shear strength. The entry and exit of critical slip surfaces at the ground surface are also applied to the model (Fig. 2). The Morgenstern-Price method, which satisfies both moment and force equilibrium equations is conducted to achieve the limit equilibrium analysis [13]. The interslice function f (x) is selected as half-sine function with a constant factor of safety distribution. The equilibrium equations used in this analysis are based on the shear mobilised at the  Fig. 2, while, the input parameters of soils and nails used in the simulation are listed in Tables 1 and 2 respectively.

2D finite difference mesh generation
A script data file using command-driven mode was developed, to generate automatically the FLAC 2D soil nailed embankment slope. This data file contains the set of FLAC 2D commands enable to discretise the geometric model shown in Fig. 1. Since the problem is plane strain model, in which lateral deformations are ignored, boundary conditions of the finite difference grid are as follows: left and right vertical sides are restrained in the horizontal direction, while at the bottom displacements are fully restricted in both directions. Furthermore, nails are simulated as one dimensional axial cable elements, which are, every nail is divided into five small parts named cable structural elements (cableSELs). On the other hand, the shotcrete surface is modeled using the structural 2D liner elements (linerSELs). Indeed, nodes are firstly created at specific locations to define the beginning and ending of these elements. Although, the nail is usually considered as a slender element, the bending effect is assumed as not important in this study. Besides, cable follows a linear elastic behaviour in tension or compression, with specified yielding limits beyond which a perfectly plastic yield maybe occurred.  On the other hand, in the soil nailing design, nails are qualified passive; they are not put in pre-tensioning at the time of their installation. It is under the effect of the soil deformations, during or after construction, and by interaction among soil and nails that the tensile nail forces have been developed. However, this interaction mechanism is very complicate. According to a simplification point of view, when the soil nail takes effect in the soil, it mainly supports the tensile stress. Indeed, the shear behaviour of nail-surrounding medium is cohesive and frictional in nature. The system is represented analogically as a spring-slider system located at the nodal point along the nail axis. In fact, the coupling spring nail element is adopted as a shown in Fig. 3. The grout annulus is assumed to behave as elastic perfectly plastic solid. Moreover, its shear behaviour during relative displacement among the nail-grout and the grout-soil interfaces as shown in Fig. 4. This behaviour is also described numerically by the grout shear stiffness K bond , grout shear strength or cohesion S bond , grout friction angle S friction and grout exposed perimeter P g [14]. It is well known that, the cohesion and friction between the soil and grout will be less than the soil effective shear strength parameters C ′ and tanϕ ′ respectively. A factor of safety (FOS) of 1.5 is therefore applied to the ultimate grout shear strength parameters in the present study [10,15]. Failure of the reinforcing system in this study is assumed to be at the grout-soil interface (Fig. 4). In such case, the shear strength of the grout parameters ( S bond and S friction ) should be evaluated at this interface as follows:

Fig. 3 Analogical model of soil nail element interaction behavior
Page 6 of 18 Derghoum and Meksaouine Geo-Engineering (2021) 12:15 where ϕ r is the friction angle of the reinforced fill [15]. The grout shear modulus G is easily given by: Usually, K bond can be directly measured from pullout laboratory tests, or it can be evaluated from a numerical estimate from the shear stress equation, describing this stress state at the grout-rock interface. In many cases, it was found that the following equation provides a reasonable estimate of K bond for the use in FLAC [16].
where G is the grout shear modulus, D is the grout hole diameter equals to 80 mm and t is the annulus thickness of the shear zone. The 2D design parameters of the soil-groutnail are listed in Table 3, while, the required input design parameters of the protective shotcrete surface for two-dimensional finite difference analysis are listed in Table 4. The corresponding numerical FLAC 2D model for nailed embankment slope is illustrated in Fig. 5.

Finite difference modelling using FLAC 3D
In this section, 3D finite difference model corresponding to the 2D model shown in Fig. 1 was generated using a specific script data file developed with the FISH language available in FLAC 3D programme. In this case, a 3D numerical model is made up of primitive zone brick elements. Nails of 12 m length are installed at 1.5 and 2 m spacing in horizontal and vertical directions respectively. A 0.15 m thick continuous elastic surface of shotcrete is also used to simulate the nail head. The embankment slope is partially resting on a soft clay of 1 m thickness (in z-direction), 14 m in depth (in x-direction) and 30 m in width Grouted-cable interaction model (in y-direction) (Fig. 6). The boundary conditions applied to the 3D model are, displacements on both vertical sides are restricted in the normal direction, while at the bottom they are fixed in all directions. Input design parameters used in 3D simulations are summarised in Table 1. However, to allow the model reaches its initial equilibrium state, only gravitational loads (self-weight) have been considered. In other words, to develop a gravitational force for this problem, a gravitational acceleration vector with a z-component of −9.81 m/s 2 has been applied in the negative z-direction. The convergence criterion,  which determines if the simulation has reached equilibrium state, is the ratio of the maximum unbalanced force to the mean gravitational body forces acting on every node (gridpoint). Although, in numerical simulation the unbalanced force never gets zero, so an acceptable lower limit maybe taken. In this study, a simulation was considered to have converged to its equilibrium state when the ratio of the maximum unbalanced force for all nodes to the average applied node forces was less than 10 −5 [17]. The protective shotcrete surface has been installed also by first creating the nodes and then creating the liner structural elements placed on the embankment slope surface as shown in Fig. 6. The shotcrete liner properties are listed in Table 5. Cable structural elements (cableSELs) are also used to simulate nails. Although, both cohesive and frictional resistances have been considered, a resisting force develops along the cable length, as response a relative motion takes place at the interface between the grout and the surrounding soil as shown in Fig. 4. Furthermore, the axial cable behaviour is linear elastic, while the grout annulus (shear zone) was assumed to behave as an elastic perfectly plastic body, so they are similar to those described for a two-dimensional model. Consequently, the shear behaviour during relative shear displacement between the nail-grout interface and the grout-soil interface can be described numerically by the grout stiffness K g (K bond in 2 D ), the grout cohesive strength C g (S bond in 2 D ), the grout friction angle ϕ g ( S friction in 2 D ) and the grout exposed perimeter P g [18]. Therefore, the 3D design parameters of the soil-groutnail are the same as those used in the 2D analysis as listed in Table 3.

Results and discussions
One of the numerical simulation advantages is the use and combination of several data. However, in order to ensure both internal and external stabilities of the nailed embankment slope, soil nailing design requires the definition of several geometric parameters. Considering the aim of this study, the results have focused only on slope stability  prediction under various soil nailing design parameters, namely nail inclination angle, nail head size and soil nail layout. Specific programmes are elaborated, to achieve perfectly these analyses and then they have been separately performed in the following sections.

Stability analysis of soil nailed embankment slope
The geometric model adopted in this section is the same as shown in Figs. 1 and 6. The effect of the nail head modes on the nailed embankment slope stability was also investigated. To simulate a fixed nail head mode numerically, a high cohesive zone is created by assigning a high grout friction ϕ g and cohesion C g to the first structural element (cableSEL). The extremity of the front end of the first part of the nail may be anchored in a specific node. In this case, there is fully compatibility between soil and nail head displacements at this node, and no stress transfer occurs around it. Moreover, a 0.15 m thick continuous elastic surface of shotcrete is used to simulate the fixed nail head mode. Unlike, when the free nail head mode is considered, the shotcrete surface is completely removed, so that there is no strong interaction between the front end of the nail and the soil close to the embankment slope surface. The choice between alternatives is important in defining the anchorages condition of nail end and the embankment slope surface, since soil mass and nail are free to move internally. Furthermore, both LEM and strength reduction methods in two (SRM2) and three (SRM3) dimensional analyses are conducted to assess the influence of this parameter on the nailed embankment slope stability. In fact, the ultimate limit state (ULS) is achieved in FLAC via the Solve Fos command. The main results of these analyses are illustrated in Fig. 7. Moreover, the shear strain rate (SSR) distribution is used to describe the slip surfaces obtained from SRM3 analysis. Although, the service limit state (SLS) refers to conditions that do not involve collapse, the ULS is progressively achieved by reducing the shear strength parameters of soils. In turn, when the free nail head is considered, the nailed embankment slope response by SRM3 at ULS, exhibits an internal face failure mode as shown in Fig. 7a. It can be seen that the tensile force in the nails increases at a constant slope from the nail head, reaches a maximum value and then decreases towards the nail head value. Besides, the line of maximum tensile forces passes behind the potential failure surfaces, but at the lower part of the embankment slope, coincides with the LEM critical sliding surface. The FOS by SRM3 is 1.53 which is slightly larger than those obtained SRM2 (1.26) , and LEM (1.17) analyses. It can be concluded that, when a lower nail-facing connection system is adopted (Fig. 7a), a facing failure mode occurs and a pullout phenomenon may then take place along the portion of soil nails located within the active zone of the embankment slope. In the other hand, when the nail head mode is used, the maximum tensile force line moves close to the embankment face (Fig. 7b). This means that, a relatively large tensile force at the nail head connection is mobilised to provide a high loading level of the embankment slope. Figure 6b shows that, when the ULS is achieved, the maximum tensile force moves slightly far from the front of embankment, in particular for the intermediate and upper rows of nails. While, it remains closer to the embankment face for the lower rows. However, failure surfaces intersect all nails indicating an internal failure mode, in which similar deeper failure surfaces are observed. It seems that they are governed by ultimate bond strength behind the slip surface (mobilised shear resistance along the soil-grout interface) and the nail length.

Overview on numerical programming process
In this section, a detailed numerical procedure to achieve a parametric study was provided by using the written programmes in the FISH programming language. Regarding the 2D finite difference modelling, the 2D numerical models were established via specific script data files that contain several command lines. For example, grid command is used to define 2D finite difference mesh limits, while generate command is then conducted to create its 2D discretised models. Propriety materials are assigned to these models via the prop command, and then boundary conditions are achieved with fix command. Again, FLAC 2D commands such as, structure node number, structure liner or cable are used to create nodes, liner and cable elements respectively at specific locations, to suit needs for the parametric study. To take account of the horizontal distance effect among the grouts (S h ) on the 2D finite difference model response, the struct prop command was used. Under which, the optional spacing keyword is specified to assign S h . Fixed nail head mode is then numerically simulated by assigning high grout friction and cohesion values to the first part of the nail via the structure chprop command. Once these programmes are completed, they can be run to check both SLS and ULS using solve and solve Fos commands respectively. In fact, 28 two-dimensional numerical models with various nail inclination angles ranged between 0 and 60° and combined with four nail lengths, namely 9, 12, 15 and 18 m, were cycled. In addition, eight other 2D models with two different soil nail layouts, considering four nail inclination angles of 0, 5, 20 and 30 • were established and cycled to design an appropriate soil nail layout. On the other hand, the 3D analysis shows more difficulties than 2D one. In which, several required subroutines were written in FISH language to perform perfectly this parametric study. As shown in Fig. 8, the programme determines the numerical model response under different nail installation angles ( α ) and lengths. It contains two main parts, in the first part two loops are used to install twenty columns ( 0 ≤ m ≤ 19) and five rows ( 5 ≤ n ≤ 9) of nails. These nails are installed at defined locations given by the beginning (x_d, y_d and z_d) and the ending (x_ed, y_ed and z_ed) coordinates. However, this programme shows that the nails are installed at 1.5 and 2 m spacing in horizontal and vertical directions respectively. It can be seen that, nail inclination angles and lengths must be defined before defining the ending x and z-coordinates to install nails with an according design. Several commands were easily added to assign nail properties or to simulate nail head mode at every step. The procedure is continued and repeated until all nails are installed. Once this step is achieved, the programme starts laying a continuous shotcrete surface. In fact, a new loop is running to define firstly the width of this surface, which is comprised between y_d and y_ed coordinates. To complete this procedure, the x and z-coordinates were updated at every step calculation. The computation process is then repeated until installation of the whole shotcrete surface. In the same way a second programme (Fig. 9) is developed to simulate the limited nail head size of 0.5 width and 0.50 height. The programme allows carrying out two limited nail heads (0 ≤ m ≤ 1) by row and three along the column (5 ≤ n ≤ 7) . The second part of this programme consists to install nails with the same order that achieved in Fig. 8. Soil nail layout with longer nail at the top can be performed following the programme shown in Fig. 10. In which, two main loops are then created with m (number of nail columns), n (number of nail rows) and len (length of nail) as variables. The process is continued and repeated until the required design was achieved. Consequently, by changing to these parameters within these main programmes, a parametric analysis can be easily carried out.

Optimization of the soil nail inclination angle
The geometric model adopted in this section is the same as shown in Figs. 1 and 6. The effect of the nail inclination angle is also investigated. The fixed nail head mode is adopted in this analysis, in which, the nail length is given for every case, while the nail inclination angle is ranged between 0 to 60 • . It was found that, when the shorter nails are used ( l nail ≤ 6m ), the influence of this parameter is negligible for both LEM and SRM analyses. Indeed, soil nailed embankment slopes with 9, 12, 15 and 18 m nail lengths, under different nail inclination angles, between 0 • and 60 • , are analysed. However, Fig. 11 summarizes the main results of these analyses. It can be seen that, two optimum nail inclination angles, of 20 • and 30 • are given by LEM and SRM analyses respectively. As shown in Fig. 11, all SRM cases give higher FOS values with an increase in ascending order. By contrast, when the LEM is used this parameter increases with a slower rate. Moreover, with the increase of the nail inclination angle, the FOS first increases and then decreases. It appears that, the FOS by SRM3 when the nail length is increased long enough ( 18 m ), becomes insensitive. Besides, whenever the nail inclination angle is then increased widely enough (> 40 • ) the FOS decreases quickly (Fig. 11d).

Effect of soil nail head size
When the soil mass is subject to deformation, the main role of the nail head is to distribute the tensile or compressive forces on the embankment slope surface. Since 3D analysis is conducted for a better consideration of soil-nail interaction in the present study, however a specific 3D model is adopted in this analysis. Indeed, for computational time reasons with a refined mesh, the width of the model as shown in Fig. 6 is reduced to 3 m , and nails of 15 m length are considered. The horizontal and vertical nail intervals are 1 and 2 m , respectively. Two different 3D model shapes of nail head modes are considered in this analysis. In the first case, a 0.15 m thick continuous elastic surface is used while in the second one the nail head size is taken as follows:0.5 m wide and 0.5 m high. In both cases, the fixed nail head process is considered. Since the computational time is very long, so for this reason only two columns and three rows of nails are considered. In this analysis, the SRM3 gives the same FOS values and critical slip surfaces for both cases (Fig. 12). This example shows that, the effect of the nail head size on the nailed embankment slope stability is not important. Although, this issue is not considered in many design guidelines, a suitable nail head size may be recommended for the stabilization of embankment slopes.

Influence of the soil nail layout
In this section, soil nail layout effect is discussed in both two and three-dimensional analyses, using two different models with fixed nail head mode. In the first case, shorter nails are installed at the top of the embankment slope, while the longer at the bottom. The nail lengths from the top to bottom are 6, 9, 12, 15 and 18 m , respectively. In the second case, the design is reversed, in which the shorter nails are introduced at the bottom and the longer at the top. Four nail inclination angles of 0, 5, 20 and 30 • are used. As shown in Table 6, FOS values given by the first design are larger than those given by the second case. It can be seen that, when the longer soil nails are in bottom, larger  Figure 13a, b show that, the critical slip surfaces by SRM3 analysis are limited at the embankment toe, while these surfaces are extended from the bottom to the top of the embankment in the second case. However, the amount of soil mass affected becomes greater (Fig. 13c, d). Therefore, case 1 is found more effective than the case 2. Another interesting issue is that the maximum tensile nail force is mobilised in the longer nail in both cases. This means that the bond stress developed in the longer nail is higher than the shorter. Since the SRM2 gives the same critical slip surfaces as those given by SRM3 analysis, only LEM slip surface with dashed line is represented in Fig. 13 for comparison. On the other hand, the shear strain ration distribution (SSR) a b c d Fig. 11 Influence of nail inclination angle on the embankment slope stability is very large at the lower part of the embankment slope. The longer nail takes loading more effectively and, hence the factor of safety is higher. Therefore, it is recommended to install the longer nails at location where failure is likely to start.

Conclusions
A numerical framework to compare soil nailed embankment slope responses from both LEM and SRM solutions, under various soil nailing design parameters is presented in this paper. Specific new programmes have been written using the FISH language built in FLAC 3D are performed in this work. Sequences of analyses are conducted and the corresponding findings can be summarized as follows: • In slope stability analysis design, the 2D software programmes are commonly used for most cases. Therefore, it was useful to compare 2D against 3D modelling results and to check how effective 2D methods are in analyzing slope design. • According to the modelling results comparison, this study has found that the FOS obtained from 2D analyses is a little conservative as compared to that of 3D. It seems that, the assumption of plane strain model without considering the horizontal limited dimension in the vertical direction of the plane is adopted in 2D analyses. Therefore, 3D analysis should be conducted to determine both FOS and potential failure surface for the embankment slope. • This study shows two failure modes, namely, a face failure will occur, when the nail head is not modeled. In turn, a pullout failure is observed when a strong nail head is considered. Although the slip surfaces by both SRM and LEM analyses usually agree well. Therefore, these failure modes were closely related to locus of the maximum tensile forces. Since the tensile force distribution, along the nail is used to provide indications on either external or internal failures, in general, the line of maximum tensile forces shown by 3D analysis does not correspond to the conventional potential failure surface as common believed. • Two different optimum nail inclination angles, namely 20 and 30° were deduced respectively from 2 and 3D analyses. Therefore, this parameter can be optimized in the soil nailing practice. • This study has found that, the effect of the nail head size on the embankment slope stability is not important. In fact, the factors of the safety obtained from this analysis are very close. • An optimum layout of the soil nail was found to be longer at the bottom and shorter at the top, which is in good agreement with some engineer's guidelines for soil nailing design during down-top construction.