the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Stability evaluation and potential failure process of rock slopes characterized by nonpersistent fractures
Wen Zhang
Jia Wang
Peihua Xu
Junqing Lou
Bo Shan
Fengyan Wang
Chen Cao
Xiaoxue Chen
Jinsheng Que
Slope failure, which causes destructive damage and fatalities, is extremely common in mountainous areas. Therefore, the stability and potential failure of slopes must be analysed accurately. For most fractured rock slopes, the complexity and random distribution of structural fractures make the aforementioned analyses considerably challenging for engineers and geologists worldwide. This study aims to solve this problem by proposing a comprehensive approach that combines the discrete fracture network (DFN) modelling technique, the synthetic rock mass (SRM) approach, and statistical analysis. Specifically, a real fractured rock slope in Laohuding Quarry in Jixian County, China, is studied to show this comprehensive approach. DFN simulation is performed to generate nonpersistent fractures in the cross section of the slope. Subsequently, the SRM approach is applied to simulate the slope model using 2D particle flow code software (PFC2D). A stability analysis is carried out based on the improved gravity increase method, emphasizing the effect of stress concentration throughout the formation of the critical slip surface. The collapse, rotation, and fragmentation of blocks and the accumulation distances are evaluated in the potential failure process of the rock slope. A total of 100 slope models generated with different DFN models are used to repeat the aforementioned analyses as a result of a high degree of variability in DFN simulation. The critical slip surface, factor of safety, and accumulation distance are selected by statistical analysis for safety assurance in slope analysis and support.
Rockslides are common geological hazards in mountainous areas. This phenomenon seriously threatens human lives and properties worldwide every year. Therefore, the analyses of rock slopes, especially their stability and failure process, are necessary for civil and mining engineering. Generally, rockslides are controlled by discontinuities (such as beddings, faults, and structural fractures). Discontinuities in rock slopes were formerly assumed in slope analysis as throughgoing and fully persistent planes. However, many researchers have suggested that discontinuities are generally not fully persistent (Eberhardt et al., 2004; Scavia, 1995; Terzaghi, 1962) and that failure surface is a combination of preexisting nonpersistent discontinuities and newly propagated cracks (Brideau et al., 2009; Einstein et al., 1983; Frayssines and Hantz, 2006; Lajtai, 1969a; Zhang et al., 2017). Nonpersistent fractures are highly complex and come in various sizes and properties, and their locations and characteristics are difficult to determine. Thus, the influence of nonpersistent fractures on slope analysis poses a great challenge (Fan et al., 2015; Wasantha et al., 2014).
Several researchers have shown that nonpersistent fractures play vital roles in the stability and failure process of fractured rock slopes (Gao et al., 2017; Huang et al., 2015; Jiang et al., 2015; Li et al., 2009; Zhang et al., 2015; Zhou et al., 2017). However, the nonpersistent fractures in these works were artificially specified by researchers as having the same sizes and orderly arranged locations. The complexity and random distribution of nonpersistent fractures in rock slopes should be considered in reflecting and modelling the stability and failure process of real rock slopes in nature.
According to a major research work on fractured rock masses, the discrete fracture network (DFN) modelling technique maximizes the use of discontinuity data from exposed surfaces and can be the best option for simulating realistic fractured rock masses (BonillaSierra et al., 2015; Chen, 2001; Elmo et al., 2013; Pine et al., 2006). By coupling the DFN technique with continuum, discontinuum, and hybrid modelling approaches, synthetic rock mass (SRM) models can be set up to investigate the mechanical properties of fractured rock masses (Elmo and Stead, 2009; Mas Ivars, 2010; Pierce et al., 2007). The SRM approach has been widely used to determine the representative elementary volume size of a fractured rock mass (Esmaieli et al., 2010), reproduce rock mass properties and behaviours (Mas Ivars et al., 2011), and simulate fracture propagation in a fractured rock mass (Zhang and Stead, 2014). SRM models have been primarily used to simulate failure and deformation of fractured rock slopes (BonillaSierra et al., 2015; Elmo et al., 2013); simulate hydraulic fracturing in naturally fractured reservoirs (Damjanac and Cundall, 2016); and estimate rock mass strength, fragmentation, and microseismicity in caving mines (Lorig et al., 2015). DFN simulation included in an SRM modelling programme presents a significant variability, which means numerously possible realizations of 2D fracture systems exist given specified input parameters (Pine et al., 2006; Zhang et al., 2020a). The results of SRM model analyses that incorporate different DFN models may vary significantly (Elmo et al., 2018; Zhang et al., 2012). Therefore, rock slope analyses based on only one DFN model can lead to erroneous results (Elmouttie and Poropat, 2014; Mas Ivars, 2006). A statistical analysis based on a large number of DFN models may reduce the aforementioned errors and provide reasonable results for rock slope analysis and support (Ferrero et al., 2016; Zhang et al., 2013, 2017).
The current study proposes a comprehensive approach that combines several wellestablished methods to conduct a stability evaluation and failure process analysis of a fractured rock slope in Tianjin city, China. First, 100 DFN models are generated on the basis of the fractures collected in the field. Second, slope models are constructed using the SRM approach. Third, the improved gravity increase method is employed to determine the stability of all slope models generated with 100 DFNs. Fourth, the potential failure processes of the fractured rock slope models are simulated. Fifth, the final critical slip surface, factor of safety, and accumulation distance are determined on the basis of the statistical analysis of 100 slope models generated with different DFN models.
2.1 Study area
A rock slope in Laohuding Quarry, which is located north of Jixian County, Tianjin, is analysed in this study (Fig. 1). The quarry was originally used for mining stromatolites, and this activity seriously damaged vegetation. After the mining area was closed, many steep slopes with internally developed discontinuities were formed. These slopes may develop rockfalls and pose potential threats to people and nearby equipment. Understanding whether rock slides will happen requires careful geological investigations.
The Laohuding Quarry area is in a lowmountain terrain, characterized by higher elevation to the north than to the south. The highest and lowest altitudes of the quarry area are 160 and 60 m, with a relative elevation of 100 m. A majority of monoclinal mountains striking south–north exist in this area. The average slopes of the mountains in the east and west of the quarry area are 25 and 30^{∘}, respectively (Fig. 1b).
This region has a continental climate with an average annual precipitation of 770.20 mm. The flood season is from June to August, and it accounts for 77.3 % of the annual precipitation. During this period, slope failures are frequently triggered. The average annual evaporation is 1867.30 mm, which is 2.42 times the amount of precipitation. These conditions indicate that the region is arid.
The lithology in this area is characterized by limestone of the middle–upper Proterozoic erathem, which exhibits a powder crystal–mud crystal structure. The limestone is moderately weathered, and the karst phenomena are scarcely visible due to low precipitation and groundwater shortage. Bedding planes, weak interlayers, faults, folds, and shear zones are not observed. This area is tectonically stable according to historic seismicity data. Nonpersistent discontinuities are randomly and widely developed in outcrops (Fig. 1c). Therefore, it is the nonpersistent discontinuities (fractures) that control the slope stability and potential failure process.
2.2 Data acquisition and DFN generation
The investigated rock slope is 20 m high and is oriented at a trend of 200^{∘} with an approximately vertical angle (Fig. 1c). Fractures in the exposed surface with trace lengths smaller than 1.5 m are widely distributed and difficult to record. Therefore, these fracture are not considered when performing fracture data collection for 2D DFN simulation. Nevertheless, the effect of these fractures on the rock mass strength is considered, which is explained in Sect. 3.1. The characteristics (such as orientation, trace length, spacing, roughness, aperture, filling, and termination) of fractures with trace lengths larger than 1.5 m were systematically surveyed using the samplingwindow method (Kulatilake and Wu, 1984). The rectangular sampling window was measured as 62 m long and 6 m high. Exactly 169 fractures were collected, and the 2D trace map is shown in Fig. 2. The mechanical properties of the different rock mass structures vary, and thus, an evaluation of the homogeneity of the fractured rock mass was conducted (Zhang et al., 2011). The result shows that this region can be considered a statistically homogeneous area. The fractures can be divided into three sets using the method proposed by Chen et al. (2005), as shown in Fig. 3. Fractures in set 1 with an average orientation of $\mathrm{39.5}{}^{\circ}/\mathrm{87.3}{}^{\circ}$ (dip direction∕dip angle) are rare. Fractures in sets 2 and 3 with an average orientation of $\mathrm{307.4}{}^{\circ}/\mathrm{44.7}{}^{\circ}$ and $\mathrm{110.2}{}^{\circ}/\mathrm{31.7}{}^{\circ}$, respectively, host the most fractures and constitute the dominant sets (Table 1). The computed Fisher constants (K) for three fracture sets are 17.1, 10.3, and 9.1, which imply the high dispersion of fracture orientations (Fisher, 1953; Priest, 1993).
The samplingwindow method features two main limits of orientation bias and trace length bias. The measured trace lengths bias occurs due to the following: the size of the sampling window is usually smaller than that of the outcrop. When the trace length is larger than the sampling window, the actual trace length cannot be totally measured; i.e., only the part of the trace length inside the sampling window can be measured. Specifically, the measured results can be divided into three types: (a) only one end of a fracture trace is measured, (b) both ends of a fracture trace are measured, and (c) no end of a fracture is measured. The mean value and probability density function (PDF) of the corrected trace lengths for each fracture set can be determined according to the methods proposed by Kulatilake and Wu (1984; Table 1). Orientation bias occurs because the fractures with small intersection angles between the fractures and exposed rock surface are more likely to be collected in the field than those fractures with large angles (Nie et al., 2020). In the present study, a 2D analysis was performed. The crosssection normal to the exposed surface was used to perform the 2D stability analysis and investigate the potential failure mechanism of the rock slope. The dip directions of fractures in set 1 are similar to that of the cross section, implying fractures acted as the surface of separation, which will not influence the results of stability analysis. Therefore, a substantial portion of the fractures of set 1 (fractures intersected by the cross section with an angle smaller than 20^{∘}) were artificially deleted prior to the stability analysis.
We derived the characteristics of the fractures (location, orientation, trace length, and density) in the cross section from those features measured across the exposed surface. Subsequently, the 2D fracture traces were generated using Monte Carlo simulation by synthesizing the aforementioned characteristics (Li et al., 2020). The specific processes are clearly described in the research of Zhang et al. (2017) and are briefly introduced as follows. The locations of the fractures were assumed to follow Poisson's distribution. We rotate the slope 20^{∘} so that the slope strikes in the NS direction and assume the fracture frequency measured along the mean normal vector direction of fracture set i is λ_{i} and the acute angle between this direction and the NS direction is η_{i}. The fracture frequency along the line parallel to the strike of the outcrop plane is λ_{i}cosη_{i} (Priest, 1993), and the crosssection plane is λ_{i}sinη_{i}. The fracture frequency of the latter is tanη_{i} times that of the former, and P_{21} (2D fracture intensity) follows this result according to the concept of the integral. The 2D orientation of a fracture is reflected by its trace gradient k, which can be expressed as sinα∕cot β (where α and β are the dip direction and dip angle of a 3D fracture, respectively). Empirical distribution was followed to generate the 3D fracture orientations based on the fracture orientation frequency, and then k in the cross section can be determined. The results of Zhang et al. (2017) revealed that the mean [E (ch^{2})] and variance [V (ch^{2})] values of the halftrace length (ch) square are constant. E (ch^{2}) and V (ch^{2}) of the collected fracture traces in the cross section can be determined according to those features on the exposed surface. The rooted numbers result in the trace lengths in the cross section. The cross section of the investigated rock slope was 20 m high, and the smallest mean dip angle of the three fracture sets was 31.7^{∘}. Thus, the length of the cross section was 30 m [$\mathrm{20}\times \mathrm{cot}\left(\mathrm{31.7}{}^{\circ}\right)$]. Finally, we used Monte Carlo simulation to merge the aforementioned parameters to generate the 2D DFN model in a 30 m × 20 m cross section.
More than one DFN can be generated on the basis of the aforementioned statistical fracture data, including the distribution type of fracture locations (Poisson's distribution), 2D fracture intensity (P_{21}), the distribution type of fracture orientation (empirical distribution), and the mean [E(ch^{2})] and variance [V(ch^{2})] values of the trace length square. For example, Fig. 4 exhibits four DFNs with the same statistical fracture data, but fracture characteristics, such as locations, dip angles, and trace lengths, are different from one another. Therefore, the rock slope analyses based on only one DFN model can lead to wrong results (Mas Ivars, 2006; Xu et al., 2014). To solve this problem, we generated and verified the validity of numerous DFN models by applying the procedure in Fig. 5. Initially, we set lines along the strike of the exposed surface, and fractures in the exposed surface can be intersected with these lines. A series of intersection points can be obtained for an individual line. Subsequently, the same procedure was performed on the generated DFNs; i.e., we intersected the DFNs using lines along the strike of the cross section and consequently generated a series of intersection points for an individual network. We compared these two sets of intersection points by using their probability density curves. When the results were identical to one another, the DFNs generated in the cross section were proved reasonable and could be selected. Finally, we selected 100 reasonable DFN models to construct different slope models and conduct stability analysis and potential failure process simulation.
The SRM approach is used to construct fractured rock slope models. In the present study, the SRM model used was a combination of two wellestablished models, namely, the bondedparticle model (BPM) in PFC2D and the DFN model (details are available in the research of Pierce et al., 2007). The SRM approach is widely used to reproduce the mechanical properties and behaviours of fractured rock masses, simulate the fracture propagation and brittle failure of fractured rock masses, and perform stability analysis of fractured rock slopes.
3.1 Parameter determination for SRM model
The SRM model in PFC2D is defined by many parameters, such as the particle contact modulus, particle normal / shear stiffness ratio, and parallelbond modulus. These parameters cannot be directly identified via laboratory and field experiments. Therefore, the parameters of the model should be predetermined according to the macroscopic characteristics of rock slopes. In particular, the input parameters of the BPM of the rock mass between the large fractures and DFN model for fractures should be ascertained and quantified separately.
Field investigation only collected the large fractures (fractures with a trace length larger than 1.5 m) which have been considered in the DFN model; thus, the effect of small fractures should be included in the BPM of the rock mass between the large fractures. The BPM parameters are calibrated against the macroproperties of the rock mass, such as the uniaxial compressive strength (σ_{c}), elastic modulus (E), Poisson's ratio, friction angle, and cohesion. These macroproperties of the rock mass are the results of comprehensive action of small fractures and intact rock material. Rock mass failure is generally shear failure, and thus the shear strength of the rock mass needs to be focused on. σ_{c}, E, and Poisson's ratio have little effect on the shear strength of the rock mass. Therefore, σ_{c}, E, and Poisson's ratio of the rock mass are assumed to be the same as those of intact rock material, which are obtained by laboratory uniaxial tests (Table 3). The friction angle and cohesion directly control the shear strength of the rock mass; thus, the friction angle and cohesion of the rock mass should be calculated by considering the effect of small fractures and intact rock material. Equation (1) is originally used to calculate the shear strength of nonpersistent fractures, which is considered to be the linear combination of rock bridges and persistent fractures (Lajtai, 1969b; Shang et al., 2018). The shear strength of nonpersistent fractures in Eq. (1) refers to that of the shear failure surface of the rock mass along the shearing direction. Therefore, Eq. (1) can be used to approximately estimate the shear strength of the rock mass. Then, the cohesion and the tangent of friction angle of the rock mass can be calculated by Eqs. (2) and (3).
where τ and σ represent the shear strength and normal stress of the rock mass between the large fractures; c_{e} and φ_{e} are the cohesion and friction angle of the rock mass between the large fractures; c_{f} and φ_{f} are the cohesion and friction angle of small fractures; c_{R} and φ_{R} are the cohesion and friction angle of intact rock material; K_{L} is the linear persistence in the shearing direction of the rock mass between the large fractures.
Laboratory uniaxial and biaxial compression tests of intact rock material obtained the values of c_{R} and φ_{R}, which are 12.25 MPa and 25^{∘}, respectively. Direct shear tests of fractures indicated the value of φ_{f} is equal to 18^{∘}. Field investigation demonstrated that no fillings existed in fractures, implying that the cohesion of small fractures is equal to zero (i.e., c_{f}=0). The linear persistence is defined as the ratio of fracture trace lengths (only the trace lengths of small fractures in the rock mass between the large fractures are considered in the present study) to the total length of the coplanar given line (Shang et al., 2018; Zhang et al., 2020b). In the present study, several lines along different directions are set in the rock mass between the large fractures, and then the linear persistence is measured. The value of the largest linear persistence measured is around 50 %. The direction of the largest linear persistence is considered the shearing direction of the rock mass between the large fractures due to the lowest shear strength in this direction. Therefore, K_{L} is equal to 0.5. Substituting the aforementioned parameter values (including c_{R}, φ_{R}, c_{f}, φ_{f}, and K_{L}) into Eqs. (2) and (3), we can obtain the cohesion and friction angle of the rock mass between the large fractures as equal to 6.125 MPa and 21.58^{∘}, respectively (Table 3).
The specimens (the heighttodiameter ratio is 2:1) for the numerical uniaxial and biaxial compression tests are set up to reproduce the macroproperties of the rock mass between the large fractures in Table 3. In PFC2D, a rock material is considered an assemblage of bonded rigid circular particles, and particle size distribution in rock models dramatically influences modelling behaviour (Mas Ivars et al., 2011). Therefore, determining the particle size of a numerical specimen remains a challenge. Generally, the ratio between the maximum and minimum radii of particles is 1.66 (Potyondy and Cundall, 2004). A small particle size is indicative of accurate simulation results. However, the number of particles of rock models increases with a decrease in particle size. A large number of particles will result in a long calculation time and low computational efficiency. In the present study, the trialanderror method was used to find a best balance of reasonable results and high computational efficiency. Specifically, we repeatedly changed the particle sizes and checked the simulation result and the time taken. Particles with radii between 0.05 and 0.083 m were finally selected to fill rock specimens.
Different macroproperties are influenced by different parameters in PFC2D. Specifically, E is mainly controlled by several parameters, including the particle contact modulus (E_{c}), particle normal / shear stiffness ratio (k_{n}∕k_{s}), parallelbond modulus (E_{b}), and parallelbond normal / shear stiffness ratio (k_{nb}∕k_{sb}); Poisson's ratio is influenced by k_{n}∕k_{s} and k_{nb}∕k_{sb}; σ_{c} is influenced by the tensile strength and cohesion of parallel bonds. Finally, the coefficient of friction is calibrated using the results of biaxial tests (Bahaaddini et al., 2013). The values of the aforementioned parameters are empirically assigned in advance, and then the numerical uniaxial and biaxial compression tests are carried out. When the macroproperties of the numerical tests correspond to the results from the laboratory tests, the parameters are considered reasonable. Otherwise, the parameters are adjusted until the rock specimens have the same macroproperties as the rock mass (Park and Song, 2009; Yang et al., 2006). The calibrated parameters of particles and bonds are listed in Table 2. The macroproperties of the rock specimens for the numerical tests and those of the rock mass between the large fractures are listed in Table 3.
In PFC2D, a smoothjoint (SJ) model is commonly adopted to simulate the large fractures in DFN. In this model, two particles that lie on the opposite sides of the intended fracture plane can overlap and interpenetrate each other instead of moving along their perimeters, thereby reproducing the real physical and mechanical properties of fractures (Bahaaddini et al., 2013). The SJ parameters are normal stiffness (${\stackrel{\mathrm{\u203e}}{k}}_{\mathrm{n}j}$), shear stiffness (${\stackrel{\mathrm{\u203e}}{k}}_{\mathrm{s}j}$), and the coefficient of friction (μ_{j}). To determine these parameters, we built a numerical specimen (with a widthtoheight ratio of 1:1) for the direct shear tests and normal deformability tests to simulate macroproperties, including shear stiffness, normal stiffness, and the friction angle. The SJ parameters can be obtained as long as the test results approximate those of the laboratory tests. Specifically, ${\stackrel{\mathrm{\u203e}}{k}}_{\mathrm{n}j}$ was obtained by the numerical normal deformability tests; ${\stackrel{\mathrm{\u203e}}{k}}_{\mathrm{s}j}$ and μ_{j} were determined by the numerical direct shear test. Table 4 lists the values of the SJ parameters and the results of the normal deformability and direct shear tests of the numerical and laboratory tests. The specific calibration procedures are complex and are thus not introduced in this paper in detail. Additional details about the calibration can be found in the works of Bahaaddini et al. (2013), Cheung et al. (2013), and Duan et al. (2016).
3.2 Model generation
The investigated rock slope with a relatively low height is located on a ground surface with a relatively flat terrain. The study area is dry; crustal movement is not obvious; and no active fault exists nearby (Sect. 2.1). Thus, the effects of ground stress, water, and earthquakes on the slope analysis were not considered in the current work.
The size of the slope section is 20 m (height) × 30 m (length). The bottom and right sides of the slope section were expanded by 10 m as the boundary section, which aims to avoid the boundary effect and does not affect the slope stability (Fan et al., 2004). Ultimately, the size of the SRM model was determined to be 30 m (height) × 40 m (length; Fig. 6). The upper boundary of the model was free. Moreover, the left, right, and bottom boundaries were assumed to be smooth rigid walls.
The particles with the same radii (0.05 to 0.083 m) as those reported in Sect. 3.1 were applied to fill the 20 m × 30 m slope section. Considering the small effect of the boundary section on the slope analysis, we filled the boundary section with particles with larger radii (0.1 to 0.15 m) to improve computational efficiency (Fig. 6). A total of 48 947 particles were generated, and the parameters presented in Table 2 were adopted for the bonds and rock particles. Notably, the values of three shear strength parameters (including the tensile strength and cohesion of parallel bonds and friction coefficient of particles) have been reduced by half according to Sect. 3.1.
Gravity was applied to the model, and the model was calculated (cycled) until the particle assemblage reached an equilibrium state (i.e., the unbalanced forces reached the required standard of 10^{−6}). Then, an embedded scripting language in the particle flow code (PFC), i.e., FISH, is used to write userdefined functions for extending the functionality or adding userdefined features in PFC. In the present study, we used the FISH functions to add the DFN into the model of the slope section by reading the location data of fractures. Subsequently, the SRM model composed of the BPM and DFN is established (Fig. 6, with the DFN in Fig. 4a as an example). The SJ parameters presented in Table 4 were adopted for the fractures. This process ignored the stress concentration at the tips of the structural fractures generated by tectonic stress, which was considered reasonable in this study since the investigated slope is characterized by the lowstress conditions and the stress concentration was intensely reduced after the longterm stability of the rock slope.
The slope was formed within one operation. Therefore, the left boundary (smooth rigid wall) was removed to simulate a onetime excavation of the slope.
4.1 Determining the factor of safety
To trigger instability, an improved gravity increase method is used in this study. The improved gravity increase method was proposed by Meng et al. (2015). This method leads to the failure of a slope in PFC2D by slowly increasing gravity acceleration and reducing the friction coefficient of particles while keeping other parameters constant. Notably, the amplitude of reduction in the friction coefficient is the same as that of the increase in gravity acceleration. In this way, the resisting force can be fixed, and therefore the factor of safety is directly reflected by the driving force. When the slope model is in a limit equilibrium state, i.e., fractures in the slope are propagated and coalesced until a throughgoing slip surface (i.e., the critical slip surface) is initially formed, the factor of safety F can be defined as the ratio of the gravity acceleration in the limit equilibrium state (g^{′}) to that in the initial state (g); i.e., $F={g}^{\prime}/g$. Taking the model in Fig. 6 as an example, we calculate the factor of safety by using the calculation procedure shown in Fig. 7. The simulation results indicate that the factor of safety of the slope model using the improved gravity increase method is 12. Additional details on the factor of safety can be found in Sect. 7 of this paper.
4.2 Initiation and propagation of fractures
When the factor of safety is determined, the critical slip surface is simultaneously obtained. To obtain an indepth understanding of the fracture propagation mechanism, the propagation process of fractures and the evolution of force chains during the formation of the critical slip surface are recorded (Fig. 8). The force magnitude is proportional to the thickness of the line in the force chain plots (lower right corner of Fig. 8), in which the blue and green colours denote compressive and tensile forces, respectively. The colour is obvious in the region where the stress concentration is strong.
Nonhomogeneous stresses are distributed throughout the slope owing to the heterogeneity of the slope model. After 2000 time steps, the compressive stress, which slowly increases from top to bottom under the action of gravity, is initially distributed throughout the slope. The tensile stress exists only in the tips of the fractures (Fig. 8a). Then, the degree of tensile stress concentration increases at the fracture end points (Fig. 8b). The fractures (black lines surrounded by a pink circle in Fig. 8) propagate from the tips of the original fractures (red lines) where the tensile stress is concentrated (Fig. 8a and b). After the initiation of fracture propagation, the tensile stress (green markings) at the tips of the original fractures dissipates (lighter green markings). A new tensile stress concentration is found at the tips of the propagated fractures, as shown in the force chain plots in Fig. 8b and c. The propagated fractures continuously expand downward to the tips of the neighbouring fractures (or rock surface) accompanied by the concentration and release of tensile stress. The orientation of the yellow arrows in Fig. 8 corresponds to the fracture propagation direction. Finally, a decrease in forces, especially in the compressive forces (lighter colour) throughout the slope, is observed; furthermore, a throughgoing surface, i.e., the critical slip surface, is formed when the propagated fractures arrive one after another at the tip of the original fractures (Fig. 8d).
According to the aforementioned results of stability analysis, the factor of safety is extremely high such that the investigated rock slope is highly stable. Nonetheless, the potential failure and the accumulation process are simulated in the present study to obtain good knowledge of the failure mechanism of fractured rock slopes and provide references for similar slope projects. To maintain coherence of analysis, we use the model in Fig. 6 as an example. Section 4 presents an analysis of the formation process of the critical slip surface of the gravityincreased model. Subsequently, another round of analysis is performed to determine the failure process after forming the critical slip surface.
Figure 9 presents the displacement field images of the slope model in different time steps. In the displacement field plots, particles are coloured according to their relative displacement magnitude. Figure 9a, which corresponds to Fig. 8d, represents the displacement field of particles when the critical slip surface is initially formed. The displacement field image indicates the occurrence of a small deformation in the failed mass above the critical slip surface; the bedrock remains stable without a distinguished gap from the failed mass. The displacements of the failed mass continuously increase, and the largest displacement of the particles (red particles in Fig. 9b) is nearly 0.2 m, which indicates the aggravation of rock failure (Fig. 9b). After 60 000 time steps, the failed mass is slowly fractured into many rock blocks along the nonpersistent fractures (Fig. 9c). The displacement and size of these blocks vary from one another. The block fragmentation is apparent near the critical slip surface because of the newly propagated fractures, whereas most of the rock blocks far from the critical slip surface remain intact (Fig. 9d). Ultimately, the failed mass is completely separated from the bedrock (Fig. 9d).
As soon as the failed mass is detached from the bedrock, a large displacement and movement of blocks occurs. To reflect the actual failure and accumulation process of rock slopes in nature, the PFC2D procedure should be manually interrupted to make corresponding adjustments (when the time steps are 80 000, which correspond to Fig. 9d). First, the particles in the bedrock are deleted to improve computational efficiency, and the critical slip surface is replaced with a rough rigid wall. The frictional properties of the rigid wall are the same as those of the propagated fractures. Then, the DFN model is removed, and the body force in the model is initialized to avoid splitting caused by the release of stress. Finally, the gravity acceleration is restored to 10.0, which corresponds to natural conditions. Figure 10 presents the computed results. In the figure, the blue sections denote the assemblage of rock particles from a macroscopic level while the red lines represent the bonds (contacts) between particles. For a convenient description, we numbered several rock blocks from 1 to 6.
In Fig. 10a (5×10^{4} time steps), except for block 2, the blocks rotate under inertia force and gravity and are separated from one another. The blocks near the slip surface are disintegrated into numerous subblocks and even crushed particles. Rock blocks 1, 3, 4, and 5 rotate anticlockwise by approximately 60, 30, 260, and 40^{∘}, respectively; whereas block 6 rotates clockwise by 75^{∘} (Fig. 10b) when the time steps are 2×10^{5}. Block 6 initially crashes to the ground, and the collision results in the bond breakage between particle blocks (red segment in Fig. 10c). Block 6 is divided into two subblocks and some crushed particles. Blocks 4, 5, and 2 successively crash to the ground one after another, but their shapes are kept intact, whereas block 3 is split owing to the collision with block 6, as shown in Fig. 10c and d. As soon as all crushed particles and blocks fall to the ground or the critical slip surface, the blocks and particles slide or roll forward as a whole and are accompanied by a fragmentation of the block edges (Fig. 10e). Figure 10f presents the final morphology of the rock model. The final deposit is composed of relatively intact rock blocks and crushed particles, and the blocks pile up above the crushed particles, presenting an inverse grading phenomenon. We also record the final accumulation distance, which is represented by the farthest reach distance of intact blocks (bonds exist in particles in these blocks). The final accumulation distance of the rock slope model shown in Fig. 6 is 28 m.
The slope analysis in Sects. 4 and 5 are based on an SRM model with one DFN model. However, as mentioned in Sect. 2, numerous DFN models can be generated because of the variability in DFN simulation, and the stability analysis of only one SRM model may lead to erroneous results. In the present study, 100 SRM models on the basis of the different DFN models (generated in Sect. 2) are built to conduct the aforementioned analysis. In particular, the critical slip surfaces and safety factors of these models are calculated using the improved gravity increase method following the calculation procedures in Fig. 7. Then, the method mentioned in Sect. 5 is employed to simulate the potential failure and accumulation process.
The factors of safety and the critical slip surfaces of the 100 slope models based on 100 different DFN models significantly vary. For example, Fig. 11 exhibits the critical slip surfaces of the SRM models based on DFN models shown in Fig. 4. According to the results, the locations and shapes of the four critical slip surfaces significantly vary, although they all extend along the nonpersistent fractures. The factors of safety of these models are 12, 15, 21, and 17.
Figure 12 illustrates the critical slip surfaces of the 100 SRM models on the basis of 100 different DFN models. Figure 13a present the factors of safety of the 100 models, ranging from 12 to 38. Therefore, the variability in the simulation results should be emphasized in the stability evaluation of fractured rock slopes. The outcome that considers numerous model calculations may lead to a rational result. In the present study, a statistical analysis method is applied to solve this problem. The final critical slip surface and factor of safety are attributed to conservative considerations. The potential critical slip surface constitutes a large failed mass (Fig. 12). The final critical slip surface covers over 90 % of the failed mass to guarantee safety in the rockslide analysis and support. Additional information on the critical slip surface can be found in the research of Zhang et al. (2017). Critical slip surfaces are supposed to have an arcshaped geometry that differs from their actual linear or broken line morphologies. Arc morphology is easily defined and is convenient for practical engineering designs. F_{s} is the final factor of safety of the rock slope when the F values of the other 90 models are greater than F_{s}. The result shows that the final factor of safety of the rock slope is 19.
Similarly, the results of the potential failure and accumulation process vary. Figure 14 presents the final accumulation states of the models shown in Fig. 11. The plots indicate that the sizes and quantities of the fractured blocks, as well as the rotating degrees, are significantly different. The final accumulation distances of the models are 28, 34, 35, and 40 m.
Finally, the different accumulation distances, with the minimum value of 15 m and the maximum value of 96 m, of the 100 SRM models based on 100 different DFN models are obtained (Fig. 13b). The potential failure and accumulation process based on one SRM model may obtain erroneous results. A statistical analysis method is also applied to solve this problem. To guarantee safety in the rock slope analysis, we attribute the final accumulation distance to conservative considerations. In particular, when the distance values of the other 90 models are lower than a certain value, then the value is the final accumulation distance (D_{a}) of the rockslide. Therefore, the final accumulation distance is 87 m.
Slope failure is a 3D stability problem. Thus, constructing a 3D SRM model for stability evaluation and failure process analysis is highly convincing and accurate. A 3D SRM model also combines rock masses in the BPM with the DFN model. However, the quantity of rock particles in a 3D slope is extremely large to conduct an effective calculation in PFC. In addition, the factor of safety obtained by 3D analysis is generally higher than that in 2D analysis. Given that many theories and technologies cannot be established and that their adoption cannot perfectly reflect rock masses at present, the accuracy of analyses may not satisfy engineering project requirements. Moreover, deriving a high factor of safety is sometimes unfavourable. Accordingly, 2D analysis, which is simple and commonly used in reality, is adopted in the present study for the stability evaluation and potential failure analysis of the investigated slope. The 2D analytical result can be regarded as a good reference, and it provides a theoretical and practical basis for future initiatives that utilize 3D analysis.
The factor of safety of the investigated slope is extremely high but reasonable. In the field investigation, weak interlayer and throughgoing discontinuities are not observed. The nonpersistent fractures are very developed, which therefore play a vitally important role in the stability of the investigated slope. The strength of the rock bridges (intact rock) is considerably higher than that of the fractures. Therefore, the obtained factor of safety is extremely high and is thus reasonable. In addition, the effects of water (rainfall) and earthquakes were ignored in the present study. However, the accuracy of the calculation result would increase if rainfall (seepage analysis) and earthquakes (kinetic analysis) are considered. This topic will be the direction of our future research.
The failure process is unlikely to occur in the investigated rock slope unless it is subjected to significant environmental changes, such as earthquakes, rainfall, unloading, or overloading. Nonetheless, the potential failure process is simulated in this study to understand the rockslide mechanism and subsequently provide a good reference for similar slope projects. For example, the size and motion of rock blocks can be utilized to predict risk degree, and the final accumulation of deposits can contribute to hazardous area division. The final arrangement of deposits (a combination of blocks and crushed particles) provides a good explanation for the inverse grading of rock avalanches reported by Cruden and Hungr (1986), Imre et al. (2010), and Wang et al. (2012).
DFN simulation presents a high degree of variability and may provide erroneous results. A statistical analysis based on numerous SRM models with different DFN models similar to those performed in the present study can reduce errors and provide conservative results for slope support. Finally, although the final results of the factor of safety, critical slip surface, and accumulation distance can guarantee safety in the rockslide analysis and support, they are not the exact values. Statistical analysis provides a new method for deriving an indepth understanding of solid earth, where specific locations and characteristics of geological materials and structures remain unknown, such as discontinuities, especially for smallscale nonpersistent fractures. Meanwhile, new theories and technologies are required to obtain precise forecasts with respect to the range values characterized by statistical methods.
The present study combines several methods, namely, DFN simulation, SRM approach, and statistical analysis, to conduct stability evaluation and potential failure process analysis of a fractured rock slope in Laohuding Quarry in Jixian County, Tianjin. The SRM technique is utilized to generate a slope model with nonpersistent fractures in the form of a DFN. The factor of safety is determined on the basis of the improved gravity increase method. The formation of a critical slip surface is also investigated. The potential failure and accumulation processes are simulated to provide a reference for similar slope projects. Numerous slope models are calculated, and the final results of the safety factor, critical slip surface, and accumulation distance are determined by statistical analysis. The major findings are summarized as follows.

The slope model with nonpersistent fractures can be effectively constructed on the basis of SRM technology. The instability of the slope model can be attained by combining the improved gravity increase method or the strength reduction method, thereby obtaining the safety factor and critical slip surface. An innovative formula to calculate the safety factor is proposed by considering stress concentration and the calculation principle of PFC2D.

Fracture propagation is closely related to stress concentration. Fractures initially propagate from the tips of the original fractures where the tensile stress is concentrated. Then, the stress is released, and a new stress concentration occurs at the tip of the propagated fractures when the fracture propagates downward to the neighbouring fractures. The critical slip surface is formed by the coalescence of preexisting fractures and newly propagated fractures.

In the initiation of failure, the failed mass is fractured into rock blocks along the preexisting fractures. Then, most blocks rotate and collapse under inertia force and gravity. Several blocks are split into subblocks owing to the collision between the blocks and the ground. The final deposit is composed of intact blocks and crushed particles, presenting inverse grading phenomena.

The critical slip surfaces, factors of safety, and accumulation distances of the slope models with different DFN models vary. Therefore, the final outcome is obtained by statistical analysis. This ensures engineering safety for rockslide analysis and support. The factor of safety (reserve) of the studied rock slope is determined to be 19. The critical slip surface is confirmed, as shown in Fig. 12. The final accumulation distance is 87 m.
The data from the research findings cannot be shared now, because these data are also part of future research.
WZ initiated this research, analysed the simulation results, and revised the manuscript. JW wrote software code and wrote this paper. PX, JL, and BS conducted field investigations and recorded field data. FW and CC dealt with the field data. XC and JQ provided guidance and collated the software code.
The authors declare that they have no conflict of interest.
This work was supported by the National Natural Science Foundation of China (grant nos. 42022053, 41877220, and 41472243), the Key Program of the National Natural Science Foundation (grant no. 41330636), the National Key Research and Development Plan (grant no. 2017YFC1501000), and the Science and Technology Project Plan of Northeast Electric Power Company (grant no. K2018N19).
This research has been supported by the National Natural Science Foundation of China (grant nos. 42022053, 41877220, and 41472243), the Key Program of the National Natural Science Foundation (grant no. 41330636), the National Key Research and Development Plan (grant no. 2017YFC150100), and the Science and Technology Project Plan of Northeast Electric Power Company (grant no. K2018N19).
This paper was edited by Filippo Catani and reviewed by two anonymous referees.
Bahaaddini, M., Sharrock, G., and Hebblewhite, B. K.: Numerical investigation of the effect of joint geometrical parameters on the mechanical properties of a nonpersistent jointed rock mass under uniaxial compression, Comput. Geotech., 49, 206–225, 2013.
BonillaSierra, V., Scholtès, L., Donzé, F. V., and Elmouttie, M. K.: Rock slope stability analysis using photogrammetric data and DFN–DEM modeling, Acta Geotech., 10, 497–511, 2015.
Brideau, M., Yan, M., and Stead, D.: The role of tectonic damage and brittle rock fracture in the development of large rock slope failures, Geomorphology, 103, 30–49, 2009.
Chen, J. P.: 3D network numerical modeling technique for random discontinuities of rock mass, Chinese Journal of Geotechnical Engineering, 23, 397–402, 2001.
Chen, J. P., Shi, B. F., and Wang, Q.: Study on the dominant orientations of random fractures of fractured rock mass, Chinese Journal of Rock Mechanics and Engineering, 24, 241–245, 2005.
Cheung, L. Y. G., O'Sullivan, C., and Coop, M. R.: Discrete element method simulations of analogue reservoir sandstones, Int. J. Rock Mech. Min., 63, 93–103, 2013.
Cruden, D. M. and Hungr, O.: The debris of the Frank slide and theories of rockslideavalanche mobility, Can. J. Earth Sci., 23, 425–432, 1986.
Damjanac, B. and Cundall, P.: Application of distinct element methods to simulation of hydraulic fracturing in naturally fractured reservoirs, Comput. Geotech., 71, 283–294, 2016.
Duan, K., Kwok, C. Y., and Pierce, M.: Discrete element method modeling of inherently anisotropic rocks under uniaxial compression loading, Int. J. Numer. Anal. Met., 40, 1150–1183, 2016.
Einstein, H. H., Veneziano, D., Baecher, G. B., and O'Reilly, K. J.: The effect of discontinuity persistence on rock slope stability, Int. J. Rock Mech. Min., 20, 227–236, 1983.
Eberhardt, D., Stead, D., and Coggan, J. S.: Numerical analysis of initiation and progressive failure in natural rock slopes – the 1991 Randa rockslide, Int. J. Rock Mech. Min., 41, 69–87, 2004.
Elmo, D. and Stead, D.: An integrated numerical modelling–Discrete Fracture Network approach applied to the characterisation of rock mass strength of naturally fractured pillars, Rock Mech. Rock Eng., 43, 3–19, 2009.
Elmo, D., Stead, D., Eberhardt, E., and Vyazmensky, A.: Applications of finite/discrete element modeling to rock engineering problems, Int. J. Geomech., 13, 565–580, 2013.
Elmo, D., Donati, D., and Stead, D.: Challenges in the characterisation of intact rock bridges in rock slopes, Eng. Geol., 245, 81–96, 2018.
Elmouttie, M. K. and Poropat, G. V.: QuasiStochastic analysis of uncertainty for modeling structurally controlled failures, Rock Mech. Rock Eng., 47, 519–534, 2014.
Esmaieli, K., Hadjigeorgiou, J., and Grenon, M.: Estimating geometrical and mechanical REV based on synthetic rock mass models at Brunswick Mine, Int. J. Rock Mech. Min., 47, 915–926, 2010.
Fan, S. C., Jiao, Y. Y., and Zhao, J.: On modeling of incident boundary for wave propagation in jointed rock masses using discrete element method, Comput. Geotech., 31, 57–66, 2004.
Fan, X., Kulatilake, P. H. S. W., and Chen, X.: Mechanical behavior of rocklike jointed blocks with multinonpersistent joints under uniaxial loading: A particle mechanics approach, Eng. Geol., 190, 17–32, 2015.
Ferrero, A. M., Migliazza, M. R., Pirulli, M., and Umili, G.: Some open issues on rockfall hazard analysis in fractured rock mass: problems and prospects, Rock Mech. Rock Eng., 49, 3615–3629, 2016.
Fisher, R.: Dispersion on a sphere, P. R. Soc. London, 217, 295–305, 1953.
Frayssines, M. and Hantz, D.: Failure mechanisms and triggering factors in calcareous cliffs of the Subalpine Ranges (French Alps), Eng. Geol., 86, 256–270, 2006.
Gao, W., Dai, S., Xiao, T., and He, T. Y.: Failure process of rock slopes with cracks based on the fracture mechanics method, Eng. Geol., 231, 190–199, 2017.
Huang, D., Cen, D. F., Ma, G. W., and Huang, R. Q.: Steppath failure of rock slopes with intermittent joints, Landslides, 12, 911–926, 2015.
Imre, B., Laue, J., and Springman, S. M.: Fractal fragmentation of rocks within sturzstroms: insight derived from physical experiments within the ETH geotechnical drum centrifuge, Granul. Matter, 12, 267–285, 2010.
Jiang, M. M., Jiang, T., Crosta, G. B., Shi, Z. M., Chen, H., and Zhang, N.: Modeling failure of jointed rock slope with two main joint sets using a novel DEM bond contact model, Eng. Geol., 193, 79–96, 2015.
Kulatilake, P. H. S. W. and Wu, T. H.: Estimation of mean trace length of discontinuities, Rock Mech. Rock Eng., 17, 215–232, 1984.
Lajtai, E. Z.: Strength of discontinuous rocks in direct shear, Géotechnique, 19, 218–233, 1969a.
Lajtai, E. Z.: Shear strength of weakness planes in rock, Int. J. Rock Mech. Min., 6, 499–515, 1969b.
Li, K. Q., Li, D. Q., and Liu, Y.: Mesoscale investigations on the effective thermal conductivity of multiphase materials using the finite element method, Int. J. Heat Mass Tran., 151, 119383, https://doi.org/10.1016/j.ijheatmasstransfer.2020.119383, 2020.
Li, L. C., Tang, C. A., Zhu, W. C., and Liang, Z. Z.: Numerical analysis of slope stability based on the gravity increase method, Comput. Geotech., 36, 1246–1258, 2009.
Lorig, L. J., Darcel, C., Damjanac, B., Pierce, M., and Billaux, D.: Application of discrete fracture networks in mining and civil geomechanics, Mining Technology, 124, 239–254, 2015.
Mas Ivars, D.: Water inflow into excavations in fractured rock – a threedimensional hydromechanical numerical study, Int. J. Rock Mech. Min., 43, 705–725, 2006.
Mas Ivars, D.: Bonded particle model for jointed rock mass, Ph.D. thesis, Royal Institute of Technology, Sweden, 2010.
Mas Ivars, D., Pierce, M. E., Darcel, C., ReyesMontes, J., Potyondy, D. O., Young, R. P., and Cundall, P. A.: The synthetic rock mass approach for jointed rock mass modelling, Int. J. Rock Mech. Min., 48, 219–244, 2011.
Meng, Y. D., Su, Q. M., Lu, W. P., and Xu, Z.: Research on improved grain flow gravity increase method, Chinese Journal of Water Resource and Power, 33, 149–151, 2015.
Nie, Z. B., Chen, J. P., Zhang, W., Tan, C., Ma, Z. F., Wang, F. Y., Zhang, Y., and Que, J. S.: A new method for threedimensional fracture network modelling for trace data collected in a large sampling window, Rock Mech. Rock Eng., 53, 1145–1161, 2020.
Park, J. W., and Song, J. J.: Numerical simulation of a shear test on a rock joint using a bondedparticle model, Int. J. Rock Mech. Min., 46, 1315–1328, 2009.
Pierce, M., Cundall, P., Potyondy, D., and Ivars, D. M.: A synthetic rock mass model for jointed rock, in: Proceedings of First Canada/United States Rock Mechanics Symposium, Vancouver, 341– 349, 2007.
Pine, R. J., Coggan, J. S., Flynn, Z. N., and Elmo, D.: The development of a new numerical modelling approach for naturally fractured rock masses, Rock Mech. Rock Eng., 39, 395–419, 2006.
Potyondy, D. O. and Cundall, P. A.: A bondedparticle model for rock, Int. J. Rock Mech. Min., 41, 1329–1364, 2004.
Priest, S. D.: Discontinuity analysis for rock engineering, London, Chapman and Hall, 1993.
Scavia, C.: A method for the study of crack propagation in rock structures, Gèotechnique, 45, 447–463, 1995.
Shang, J., West, L. J., Hencher, S. R., and Zhao, Z.: Geological discontinuity persistence: Implications and quantification, Eng. Geol., 241, 41–54, 2018.
Terzaghi, K.: Stability of steep slopes on hard unweathered rock, Géotechnique, 12, 251–270, 1962.
Wang, Y. F., Cheng, Q. G., and Zhu, Q.: Inverse grading analysis of deposit from rock avalanches triggered by Wenchuan earthquake, Chinese Journal of Rock Mechanics and Engineering, 1089–1106, 2012.
Wasantha, P. L. P., Ranjith, P. G., Xu, T., Zhao, J., and Yan, Y. L.: A new parameter to describe the persistency of nonpersistent joints, Eng. Geol., 181, 71–77, 2014.
Xu, Q., Chen, J. Y., Li, J., and Yue, H. Y.: A genetic algorithm for locating the multiscale critical slip surface in jointed rock mass slopes, Math. Probl. Eng., 2014, 543081, https://doi.org/10.1155/2014/543081, 2014.
Yang, B., Jiao, Y., and Lei, S.: A study on the effects of microparameters on macroproperties for specimens created by bonded particles, Eng. Computation., 23, 607–631, 2006.
Zhang, K., Cao, P., Meng, J. J., Li, K. H., and Fan, W. C.: Modeling the progressive failure of jointed rock slope using fracture mechanics and the strength reduction method, Rock Mech. Rock Eng., 48, 771–785, 2015.
Zhang, W., Chen, J. P., Liu, C. H., Ma, J. Q., and Zhang, C.: Application of chisquare test in structural homogeneity dividing of fractured rock mass, Chinese Journal of Geotechnical Engineering, 33, 1440–1446, 2011.
Zhang, W., Chen, J. P., Liu, C., Huang, R., Li. M., and Zhang, Y.: Determination of geometrical and structural representative volume elements at the Baihetan Dam Site, Rock Mech. Rock Eng., 45, 409–419, 2012.
Zhang, W., Chen, J. P., Zhang, W., Lv, Y., Ma, Y. F., and Xiong, H.: Determination of critical slip surface of fractured rock slopes based on fracture orientation data, Sci. China. Technol. Sc., 56, 1248–1256, 2013.
Zhang, W., Zhao, Q. H., Chen, J. P., Huang, R. Q., and Yuan, X. Q.: Determining the critical slip surface of a fractured rock slope considering preexisting fractures and statistical methodology, Landslides, 14, 1253–1263, 2017.
Zhang, W., Fu, R., Tan, C., Ma, Z. F., Zhang, Y., Song, S. Y., Xu, P. H., Wang, S. N., and Zhao, Y. P.: Twodimensional discrepancies in fracture geometric factors and connectivity between fieldcollected and stochastically modeled DFNs: A case study of sluice foundation rock mass in Datengxia, China, Rock Mech. Rock Eng., 53, 2399–2417, 2020a.
Zhang, W., Lan, Z. G., Ma, Z. F., Tan, C., Que, J. S., Wang, F. Y., and Cao, C.: Determination of statistical discontinuity persistence for a rock mass characterized by nonpersistent fractures, Int. J. Rock Mech. Min., 126, 104177, https://doi.org/10.1016/j.ijrmms.2019.104177, 2020b.
Zhang, Y. and Stead, D.: Modelling 3D crack propagation in hard rock pillars using a synthetic rock mass approach, Int. J. Rock Mech. Min., 72, 199–213, 2014.
Zhou, Y., Chai, J. F., and Han, G.: Meso numerical study on the failure mechanism on rock slope with bedding intermittent joints, Engineering Review, 37, 272–280, 2017.
 Abstract
 Introduction
 Study area and data acquisition
 SRM model for rock slope analysis
 Stability analysis
 Potential failure and accumulation process
 Statistical analysis
 Discussion
 Conclusions
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Study area and data acquisition
 SRM model for rock slope analysis
 Stability analysis
 Potential failure and accumulation process
 Statistical analysis
 Discussion
 Conclusions
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References