#### Abstract

The mudded weak interlayer is a geotechnical sandwich material exhibiting strain softening behavior, which plays an important part in the slope stability. The present work primarily focuses on the shear strength of the mudded weak interlayer in rock slopes. To determine the peak and residual shear strengths of the mudded weak interlayers, the particle flow code (PFC) is used to simulate the failure behavior during the direct shear tests. Laboratory investigations including uniaxial compression test, SEM, and 3D deformation measurement are employed to calibrate the essential micro parameters of the mudded weak interlayer during the simulation process in PFC. The numerical model is built based on these parameters and both the peak and residual shear strengths can be predicted from the model. The prediction results show that the peak and residual internal friction angle are 19.36° and 14.61°, while the peak and residual cohesion are 22.33 kPa and 2.73 kPa, respectively. Moreover, to validate the obtained peak and residual strengths, the results are compared with literature data. The peak and residual shear strengths of the mudded weak interlayer can serve as an important benchmark to evaluate the stability of side slopes and provide guiding suggestions for their reinforcement.

#### 1. Introduction

The mudded weak interlayers are thin and weak rock strata in slopes, which are formed by tectonism and weathering. Because the shear strength of these weak interlayers is lower than that of the other layers in slopes, they are considered to be the most dangerous potential slip surfaces [1–4]. Thus, the methods for the evaluation of the shear strength of the mudded weak interlayers have been the focus of research [5]. The dynamic change of the shear strength of the mudded interlayers from the peak to the residual was often neglected in the traditional analysis of the rock slope stability [6]. However, with the increase of shear deformation, strain softening will occur in the mudded weak interlayers [7–9]. In view of the importance of studying the shear behavior of the mudded weak interlayers in slope stability evaluation, it seems essential to develop a method that reasonably predicts the peak and residual strength of the mudded weak interlayers.

Extensive attention has been paid to the estimation of the shear properties of the mudded weak interlayer in recent years [10, 11]. More specifically, a thorough survey shows that the existing studies mainly concentrate on two aspects, the general shear properties and the strain softening characteristics.

Firstly, soil tests, including the large ring shear tests and uniaxial and triaxial compression tests, were employed to investigate the general shear properties of the mudded weak interlayers. Li et al. [12] investigated the mechanical behavior of bedded rock salt containing an inclined interlayer using experimental techniques including uniaxial and triaxial compression tests. The research demonstrated that the inclined interlayer has an effect on the deformation and fracture of rock salt [12]. Xu et al. [13] developed a simple shear strength model for the interlayer shear weakness zone based on the available experimental data, which is capable of describing the shear behavior of both the interlayer soil and the soil/rock interface. Li et al. [14] studied the shear properties of slip zone soils from three giant landslides in the Three George Project area using the large ring shear tests at different shear rates (0.1, 1, and 10 mm/s). According to the test results, the relations between the Atterberg limits, particle size distribution, and shear properties were investigated. Chen et al. [15] studied the physical and mechanical properties of the mudded weak interlayer in a natural slip zone of red beds in southwestern China. It was found that the mudded weak interlayer of the progressive landslide is composed of fine-grained soils with a large amount of clay particles and the mudded weak interlayers are overconsolidated.

The study on the strain softening characteristics seems more popular than that of the general shear properties of the mudded weak interlayers. Terzaghi et al. [17] first presented the concept of strain softening and concluded that the fissured clay has a higher shear strength when the shear strain value is at a low level. Based on direct shear tests of overconsolidated soils, Skempton [18, 19] found that the shear strength of the overconsolidated soils will decrease when the shear strain exceeds a certain value. Using an extended Mohr–Coulomb constitutive model to represent the strain softening behavior of slope material, Mohammadi and Taiebat [7] evaluated the postfailure deformation of slopes and embankments with a numerical method based on the updated Lagrangian formulation. The results showed that the postfailure deformation of the slopes is a function of the strength reduction rate and the stiffness of the slope material. Conte et al. [8] also presented a numerical approach to analyze the stability of slopes in soils with strain softening behavior. The strain softening behavior of the soil was simulated by reducing the strength parameters with the increasing deviatoric plastic strain. Chen et al. [15] reported that the mudded weak interlayer exhibits a strain softening behavior, based on the shear stress–strain curves and that the residual shear strength is closely related to the water content. The microstructure of the mudded weak interlayer was also investigated using polarized microscopy and scanning electron microscope techniques. Wen et al. [20] investigated the relation between the residual strength and the index properties (particle size distribution, Atterberg limits, and mineral composition) of slip zone soils for 170 landslides in the Three George Project area. The relation can be used to evaluate the residual shear strength of the slip zone soils. Using laboratory triaxial tests, Indraratna et al. [21] found that the peak shear strength of the rock joints with compacted infill (the interlayer) increased with the decrease of the degree of saturation from 85% to 35%. Based on the laboratory observations, they developed an empirical model to describe the shear strength of the rock joints with the weak interlayer. Papaliangas et al. [22] concluded from a series of direct shear tests on sandstone discontinuity models with mean roughness amplitude of 7 mm and filled with weak interlayers of varying thickness that the peak shear strength will decrease with the increase of the thickness of the weak interlayers; however, the residual shear strength will decrease less markedly.

Despite the fact that the existing studies [7, 8, 12–15, 17–23] have presented available experimental techniques to investigate the shear behavior of the mudded weak interlayers and the softening strain characteristics, certain unavoidable problems still challenge researchers during the determination of the shear strength. For example, a small-scale sample may differ from the actual rock slope deformation behavior due to the scaling effect. The best method is the full-scale study; unfortunately, this is costly and difficult to execute. Therefore, the traditional laboratory tests and in situ tests have limitations. In consideration of the significance of the determination methods for the shear strength of the mudded weak interlayers and the limitations of the existing testing methods, it is necessary to develop a new approach to study the shear behavior and the strain softening characteristics of the mudded weak interlayers. Although the particle flow code (PFC) is widely used for geotechnical engineering materials [24–28], the predictive models for the shear strength of the mudded weak interlayers based on the PFC have not yet been reported in the literature, to the best of our knowledge.

DEM, specially applied to solve the discontinuous medium-related mechanical problem, is suitable for the simulation of sandwich material in slopes. On one hand, the sandwich materials in slopes, although thin and weak, are made of large numbers of weathered rock particles from the perspective of mesostructure. On the other hand, the mechanical simulation based on the continuum mechanics theory cannot connect the macroscopic properties of materials with their microstructure characteristics. However, the DEM simulation, such as PFC, is able to not only reflect the interactions among the particles, but also obtain the macroscopic mechanical properties. In this paper, we aim to determine the peak and residual shear strengths of the natural mudded weak interlayers. For this purpose, a numerical model for the mudded weak interlayers is established using the PFC in two dimensions (PFC 2D), which is a DEM numerical method. To match the numerical model better with the real conditions of the mudded weak interlayers, the mesomechanical parameters were calibrated using laboratory tests. Specifically, due to the limitations of the direct shear and triaxial compression tests (such as the dimensional limit), the uniaxial compression tests and deformation monitoring using digital speckle measurement techniques (ARAMIS) were applied to adjust the numerical model parameters to the optimal values. Thus, the strain softening characteristics of the mudded weak interlayers can be determined through a series of numerical shear tests and the shear strength (including the peak and residual strength) can be obtained from the test results. Finally, the determination method for the shear strength of the mudded weak interlayer presented in this paper is validated by comparison of the prediction results with the existing literature values.

#### 2. Laboratory Investigation

To study the shear behavior of the mudded weak interlayers and reasonably predict the shear strength, the PFC 2D was employed to establish the numerical forecast model based on which both the shear failure mechanism and shear strength can be obtained. The mesomechanical parameters of the PFC models cannot be directly measured with experimental methods; instead, uniaxial compression tests have to be carried out on the mudded weak interlayer specimens to calibrate the microparameters of the particles. Deformation monitoring of the specimen surfaces using digital speckle measurement techniques (ARAMIS) provides additional information, which is useful for the calibration.

##### 2.1. Study Area and Specimen Preparation

Figure 1 is a global image of the mudded weak interlayers (a) and the collected mudded weak interlayer specimens (b). All of the specimens were obtained from the right side of the slope of tunnel number 2 of the Longquan mountain, Southwest China, where red mudstone is widely developed. Both the overlying rock strata and the underlying bed of the mudded weak interlayers are red mudstones (Figure 1(a)). To keep the structure and the water content of the specimens undisturbed, all the mudded weak interlayer specimens were collected by cutting ring samples, which were sealed with wax (Figure 1(b)). Each cylindrical specimen has a diameter of 61.8 mm and a height of 20 mm.

**(a)**

**(b)**

##### 2.2. Test Procedure

The uniaxial compression tests for the mudded weak interlayers were carried out using a triaxial apparatus, applying only the axial loads. As mentioned above, the mudded weak interlayer is a kind of thin and weak rock strata, so it is not that likely to obtain a sample that meets the dimension standard (2 : 1 ratio). The authors employed the digital speckle ARAMIS to make an aided analysis of the UCS results. Specifically, the three-dimensional optical measurement system (ARAMIS 3D) was applied to monitor the strain distribution along the surface of the specimen under the axial load. ARAMIS is capable of recording the deformation of the specimens in the uncontact condition with a high precision of 0.001 mm. ARAMIS should be calibrated before the experiment to ensure the precision of the results (Figure 2).

**(a)**

**(b)**

Subsequently, the surface of the specimens should be pretreated with digital speckles. The CCD of ARAMIS can obtain the images of the specimen surface with the change of the axial load and transform them to digital signals. The grey scale information of the randomly distributed speckle field on the surface of the object is tackled using the digital image correlation method and the image is divided into a grey-level matrix. The deformation of an arbitrary point on the specimen surface can be measured by investigating the movement of the corresponding subarea. Figure 3 shows the deformation of the subarea on the surface of the specimen.

The point is used as the center of the selected subarea of the reference image (undeformed) and (, ) is the center of the subarea of the target image (deformed; Figure 3). Based on the assumption that the point () is another point in the selected subarea of the reference image, the new coordinate of the point (, ) in the target image can be described as follows:where and denote the displacements of the point in the directions of and .

As mentioned above, all of the monitored specimens must be specially processed, as shown in Figure 4. Firstly, the surfaces of the specimens are sprayed with white paint and are treated as the background. After air-drying of the white paint, the speckles were uniformly sprayed on the specimen surfaces with black paint.

**(a)**

**(b)**

Three groups of mudded weak interlayer specimens with special speckles on the surfaces were tested using the hydraulic loading system REY-8000. The maximum axial load of this loading system can reach 100 KN and the maximum displacement is 530 mm. During the entire test, the axial load and displacement were recorded using the axial load cell and the displacement transducer, respectively. Meanwhile, the deformation of the specimen surfaces was recorded with ARAMIS 3D. All of the testing results were used to calibrate the micromechanical parameters of the mudded weak interlayer model in the PFC2D. The loading framework in the displacement-control mode and the loading rate are 5 mm/min, which is higher than the common loading rates of the UCS test. However, the ARAMIS 3D is able to capture the pictures of the specimens with a high frequency. The smooth plastic slices were added to the upper and lower surface of the specimens during the tests, as shown in Figure 5, to eliminate the influence of friction.

**(a)**

**(b)**

##### 2.3. Testing Results and Discussion

###### 2.3.1. Typical Stress–Strain Curves from Uniaxial Compression Tests

The failure criterion of the mudded weak interlayers is extremely important for the determination of the maximum vertical stress, which can be applied to the calibration of micromechanical parameters of the PFC models. In stress space, the ideal plastic yield surfaces are supposedly invariable, and the initial yield surfaces are assumed the same as the subsequent yield surfaces [29, 30]. Therefore, no criteria are available in stress space to determine the failure of geotechnical materials. However, the criteria in strain space can overcome these deficiencies and the permissible deformation can serve as an intuitionistic failure criterion for the mudded weak interlayers, which can be described as follows:where denotes the maximum permissible strain of the mudded weak interlayer. The value of the maximum permissible strain can be determined by experiments. As presented in [31], the normally consolidated clay was investigated using drained triaxial compression tests with the confining pressures of 70 kPa, 200 kPa, and 700 kPa. The resulting stress–strain curves show that the axial deviatoric stress continues to increase until the strain exceeds 20%, which can be chosen as the maximum permissible strain .

In this paper, the maximum permissible strain of the mudded weak interlayer can be determined using a similar approach. Because of the limited amount of mudded weak interlayer specimens available for the uniaxial compression tests, three typical specimens were tested. Figure 6 shows the typical stress–strain curves during the uniaxial compression. The deformation process of the mudded weak interlayer upon increase of the axial load can be divided into two stages (Figure 6). In the first stage, which is also known as the elastic stage, the axial strain of the specimens increases linearly with the increasing load. The slope coefficients of the curves mainly depend on the stiffness of the specimens in this stage and the mudded weak interlayer is not destroyed yet. The second stage is known as the plastic stage. The axial load stays constant when the strain exceeds 10%, while the strain keeps increasing (Figure 6). Therefore, the strain value of 10% can be chosen as the maximum permissible strain of the mudded weak interlayer and the axial load when the strain reaches can be treated as the maximum vertical stress (peak stress; Table 1).

The average maximum vertical stress of the three mudded weak interlayer specimens has been obtained as follows:

###### 2.3.2. Major Strain Monitoring Results

Because the loading rate is 5 mm/min and the height of the specimen is 20 mm, the loading time is less than 30 s. According to the actual demands of deformation monitoring, the photograph frequency was set to 5 times per second. To record the dynamic response process of the surface strain changing with the vertical stress, the entire deformation monitoring procedure was divided into 7 loading stages. Due to the little strain of the specimen at low loads, which is inconspicuous, the loading time of first stage is 6 seconds. The loading time for each of the second to seventh stages reduces to 3 seconds with the increase of the vertical stress and strain. Based on the mesh generation and unit analysis of the monitored images, the evolutionary trends of the major strain for the mudded weak interlayer specimen under uniaxial compression are shown in Figure 7.

Figure 7 shows that local deformation occurs with the increase of the vertical load. It is worth noting that the surface deformation of the specimen is inhomogeneous. The unevenly distributed deformation of the specimen surface is mainly caused by the anisotropy of the inner structure of the mudded weak interlayers. As the load increases further, the plastic zone develops continuously and the shear band begins to occur in the locally weak area. The shear band keeps developing along the longitudinal direction until it crosses the whole specimen.

More specifically, the deformation evolution of featured sections and points was obtained from the data collecting and analyzing system of ARAMIS. As shown in Figure 8(a), three sections were selected for the investigation of the local deformation characteristics of the specimens. The section monitoring results are shown in Figure 8(b). Figure 8(b) shows that all the units of Section of the last loading stage are at a high level of strain, with the major strain > 10%. As a result, Section can be treated as the failure surface. In contrast, the major strain of the units in Section remains at a low level and the maximum strain is <7.5%. Finally, the major strain for the units with the same vertical position but different horizontal position in Section is also shown in Figure 8(b).

**(a)**

**(b)**

Similarly, four feature points were set on the surface of the specimen to monitor the dynamic changing process of the major strain. Figure 9(a) shows the position of the feature points among which the points 0 to 2 are in the shear band of the mudded weak interlayer and point 3 is in the center of the specimen. The values of the major strain for the points 0 to 2 keep growing linearly during the loading process (Figure 9(b)). However, that of point 3 hovers at zero, which means that the major strain remains stable.

**(a)**

**(b)**

###### 2.3.3. Prediction of the Cohesion and the Internal Friction Angle of the Specimens

As known from classical soil mechanics, the failure envelop cannot be obtained because only one Mohr circle of stress can be plotted based on the uniaxial compression testing results. According to the Mohr–Coulomb intensity limit equilibrium of soil, the cohesion and the internal friction angle can be determined only if the break angle of the specimen can be measured. For the unconfined compression test, the limit equilibrium condition can be described as follows:where equals the maximum vertical stress and the break angle can be expressed as 45° + . Thus, the cohesion and internal friction angle can be obtained as follows:

There are two notable breaking sections and the break angles can be automatically measured by ARAMIS (Figure 10). The final break angle is determined using the method of averaging based on which the cohesion and internal friction angle can be obtained as follows:

It is notable that the cohesion and the internal friction angle are obtained from the Mohr–Coulomb intensity limit equilibrium, which can be regarded as the peak shear strength. However, the residual shear strength of the mudded weak interlayers is still unknown. The test results, including the cohesion and the internal friction angle, are used to calibrate the microparameters of the bonded particles in PFC2D. Based on the calibrated model for the mudded weak interlayers, the strain softening characteristics can be simulated while both the peak and the residual shear indexes can be predicted.

#### 3. Numerical Model for the Mudded Weak Interlayers Based on the Test Results

##### 3.1. Numerical Model Set-Up

###### 3.1.1. Bond Style Selection

Both the contact bond (CB) and parallel bond (PB) are provided in PFC [32]. The PB model can be understood as a circular (2D) or rectangular (3d) cross-section of the contact plane, which transmits force and momentum among the particles. However, the CB model approximates the physical behavior of cement similar to a spring with both the normal and tangential bond strengths, which are only active at the contact point and can only transmit force [32]. The damage of the aggregated particles is reflected in the breakage of the bonds. Compared to the PB model, the CB model has been more widely used because of its reduced number of microparameters [10]. Therefore, the CB model is used in the current study. To generate a numerical model for the simulation of the mudded weak interlayers using PFC 2D, the micromechanical parameters are calibrated against the test results as described above.

###### 3.1.2. Basic Geometric Parameter Design

According to the actual size of the specimens, the numerical sample has a height of 20 mm and a width of 61.8 mm (Figure 11). During the direct shear test in PFC 2D, eight walls are defined. Subsequently, 2238 particles are generated inside the defined walls after setting the maximum and minimum diameters ( mm, mm). Undoubtedly, as DEM numerical software, PFC2D has its own limitations; uncertainty factors may result in different test results, such as that of the particle size and damping set. However, it has been proven to be a useful and reliable tool to study geotechnical materials. This study attempts to optimize the numerical model of the mudded weak interlayers in combination with the laboratory test results and other techniques.

###### 3.1.3. Determination of the 2D-Porosity

It is important to note that the PFC model of the mudded weak interlayer is two-dimensional. However, the traditional porosity tests, such as the pycnometer method, only test the 3D porosity. In this paper, the 2D-porosity of the specimen is obtained from the image segmentation of the mudded weak interlayer section micrograph. It is worth mentioning that all the section micrographs of the mudded weak interlayers are captured using the electronic scanning microscope. To eliminate the image noise and enhance the edges, spatial domain enhancement [33] was adopted to address the original section micrographs. This method can reduce the noise of the image effectively while it enhances the contrast of the image (Figure 12). Finally, threshold segmentation is applied to process the section micrographs, which is typically used to locate objects and boundaries (lines or curves). Thus, the solids and pores are separated into two independent parts, which are shown in Figure 13. The porosity can be determined using MATLAB. Mathematically, the 2D-porosity of the whole section micrograph is the statistical average of the local section micrographs in the section micrograph. According to the test results, the porosity of the numerical model is set to 0.279.

**(a)**

**(b)**

**(a)**

**(b)**

###### 3.1.4. Setting of Wall Stiffness

In PFC, the force cannot directly act on the walls. Instead, the servomechanism is applied to control the stress of the walls. Both the loading process and the model constraints are accomplished through the walls. Therefore, the setting of the wall stiffness is of great importance to the accuracy of the simulation results. The overlarge wall stiffness would increase the time steps of the initial equilibrium state, while the walls would be pierced by the balls if the wall stiffness is too small. The stiffness of the servowalls is therefore set to one-tenth of the ball stiffness to simulate the flexible boundary of the servowalls. However, the stiffness of the loading walls is ten times that of the ball stiffness.

##### 3.2. Calibration of the Contact-Bonded Model

In the contact-bonded model of the mudded weak interlayer in PFC 2D, the unknown parameters mainly include the particle stiffness (normal contact stiffness and shear contact stiffness ), the friction coefficient, the normal and shear strengths of the contact bond (NB and SB), and the loading rate . To calibrate the micromechanical parameters of the numerical model more efficiently, the sensitivity between the macro- and microparameters needs to be analyzed before the calibration process. For this purpose, the influence of the microparameters on the macroproperties can be investigated by changing the values of the mesoscopic structural and mechanical parameters. To make full use of the experimental test results, both the direct shear and uniaxial testing models are established to perform the sensitivity analysis. The macroscopic response mechanism upon change of the microscopic parameters is deduced from the numerical test results (Table 2).

The modulus of elasticity mainly depends on the normal contact stiffness and the shear contact stiffness of particles, as well as their ratio . The macrostrength parameters are closely related to the normal and shear bond strengths (NB and SB) and are affected by the friction coefficient among the particles. Meanwhile, the friction coefficient has an important effect on the residual strength of the numerical model. In addition, the porosity , vertical stress , and the loading rate have a disproportionate level of impact on the macroproperties, which can be seen in Table 2. All the microscopic properties of the particles and bonds are calibrated by trial and error. More specifically, the microparameters are adjusted iteratively to match the numerical testing results to the experimental results. Table 3 shows the calibrated values of the main microparameters at different vertical stress.

#### 4. Strain Softening Characteristics and Numerical Prediction of Peak/Residual Shear Strength

As discussed in Section 3.2, all of the essential microscopic parameters are obtained through the calibration based on the laboratory test results. As a matter of fact, the calibrated model of the mudded weak interlayer is available for the numerical tests. Direct shear tests under constant vertical pressures of 50 kPa, 100 kPa, 150 kPa, 200 kPa, 250 kPa, and 300 kPa are simulated using PFC 2D and the results are presented in Figure 14.

Figure 14 shows that the shear stresses increase almost linearly with the shear strains before the stress peaks in the first stage, which is also called the linear elastic stage. The stresses begin to slip down after that, while the shear strains keep rising until they reach a relatively stable stage. Actually, because the physical properties of the mudded interlayers are complex, the functional curves are converging slowly. Thus, the least-squares fitting method is used to obtain an approximate stable value (Figure 15). Figure 15 shows the functional relations between the shear stresses and shear strains of the data plots after the peak point from binomial fitting. The shear stresses gradually approach a steady value, which is the residual shear strength. Both Figures 14 and 15 show that the peak shear strength (PS) and the residual shear strength (RS) of the mudded weak interlayers keep an upward tendency with the increase of the vertical pressure. At the same time, the corresponding strains also increase at different levels. It is notable that the ratio between PS and RS stays in the range from 0.5 to 0.7.

According to the direct shear test results under different vertical pressures, PS and RS are plotted as functions of the normal stresses in Figure 16. Least-squares fitting is applied to obtain the functional relations between , , and the vertical pressure :

It can be determined from (7) that

As mentioned in Section 2.3, the friction coefficient and the cohesion of the mudded weak interlayers estimated by the uniaxial compression test and deformation monitoring is and . Compared with the numerical predicted peak friction and the cohesion, the similarity index is 96.8% and 90%, respectively.

#### 5. Verification of the Present Prediction Method

The PS and RS of the mudded weak interlayer can also be obtained by repeated direct shear tests [6]. Hu and Zheng [6, 16] surveyed the PS and RS indexes of the common mudded weak interlayers, which appear frequently in mudstone, carbonatite, and carbonaceous shale. In this paper, the mother rock of the mudded weak interlayer is red mudstone.

Both the numerically predicted shear strength indexes of the mudded weak interlayer and the values obtained by statistics are presented in Table 4. The peak and friction coefficients determined by the present method are in the range of the statistics. As for the cohesion, the numerically predicted result is a little bigger than the regular values. In general, the predicted strength parameters are in good agreement with the statistic values from the literature. These facts can be used to validate the feasibility of the present method in predicting the shear strength parameters of the mudded weak interlayers.

#### 6. Conclusions

In this study, we presented a method to study the strain softening characteristics of the mudded weak interlayers and predict their shear strength indexes. For this purpose, soil experiments, deformation monitoring techniques (ARAMIS), and DEM numerical simulation are employed. To simulate the real state of the mudded weak interlayers, the microscopic parameters of the numerical model are calibrated based on the soil test results. Specifically, the deformation monitoring during the uniaxial compression tests acts as an important benchmark in the calibration process. Furthermore, the calibrated model is used to perform several numerical direct shear tests at different vertical pressures. The shear strength indexes of the mudded weak interlayers are determined using the numerical calculation results. From the comparison of these predicted shear strength indexes with the statistic values, the validity of the present method can be verified. The relatively accurate PS and RS indexes of the mudded weak interlayers can serve as important technical tools to evaluate the side slope stability.

#### Conflicts of Interest

The authors declare that they have no financial and personal relationships with other people or organizations that can inappropriately influence their work; there is no professional or other personal interest of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in, or the review of, this manuscript.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (no. 51574201) and the State Key Laboratory of Geohazard Prevention and Geoenvironment Protection (no. SKLGP2015K006).