The assessment of the strength parameters of geological formations in regional scale which encounters thousands of slopes is a complicated approach and time consuming and needs huge field work. This issue is an important research topic concerning the regional seismic-landslide susceptibility analysis or hazard zonation. An empirical regression model was presented to estimate the Geological Strength Index (GSI) with an implication on geological quadrangle of Gorgan region at Alborz mountains range (north of Iran). Two main sets of data were applied in this study: (1) geomorphological data including the slope height, aspect, and distance from faults and distance from thrusts and (2) the physical and mechanical properties of rocks including the unit weight, uniaxial compressive strength (σci), and the petrographic constant (mί) of intact rock. The first group was extracted from a 1 : 100,000 digital geologic map and 10 m digital elevation model (DEM) and the second group was obtained from the Hoek–Brown failure criterion recommended tables. Linear regression equations were generated applying data collected from 294 studied stations using SPSS software. The regression equation predicted GSI in terms of (1) the distance from faults, (2) the distance from thrusts, and (3) the uniaxial compressive strength (σci). The equation had an R2 value of 0.739 and thus fit well to the data. The new method in its present state was recommended for the estimation of the GSI values in regional scale conditions for the assessment of landslide susceptibility and hazard mapping or post events landslide occurrence prediction in the case of probable big earthquakes in Alborz area that is required for emergency responses. The results indicated that the estimation error was about ±30 while the average error was within +5 and −5 and average error percentage was about 3%.

1. Introduction

Since the last decade, using Newmark’s displacement method [16] accompanied by the Geographic Information System (GIS) has become a useful deterministic approach to prepare the seismic-landslide hazard map. Although many general studies have been conducted on the identification and description of landslides in general, the application of GIS based on Newmark’s method is relatively new in Iran, and these types of studies have been just carried out in recent years [7]. Mahdavifar [8] generated a fully automatic version of a GIS-based system based on a simplified Newmark’s displacement method which could provide a seismic-landslide hazard zonation map immediately after the occurrence of an earthquake, [2, 4]. This study was conducted in Alborz and central areas of Iran. In this method, hazard zonation is performed by calculating safety factor, using slopes steepness indices, material parameters (shear strength), and the characteristics of the expected earthquake. Accordingly, in this method, the shear strength of the geological units is a fundamental issue.

In recent years, the Geological Strength Index (GSI) classification system, proposed by Hoek [9] and Hoek et al. [10], has been considered as the most acceptable empirical method to estimate the rock mass strength and deformation parameters, and it can be said that there are no other suitable alternatives for it [11, 12]. The GSI is a technique to estimate the rock mass strength in different geological conditions, using some standard charts, site observation, and rock mass description. The rock mass properties can be determined considering the degree of crushing and conditions of discontinuities surfaces, indicated by roughness and alteration. Combining these two parameters provided a principal basis to describe rock mass types with diversified rock structure ranging from very tightly interlocked strong rock fragments to heavily crushed rock masses. Based on the rock mass description, the GSI value could be estimated by the contours to reach a value of 0–100, representing the overall geotechnical quality of rock masses. Since the last decade, the GSI approach has been modified by many researchers [1217]. Bieniawski [18] and Hoek and Brown [10] suggested the relationship between GSI and rock mass quality index Q and rock mass rating (RMR), respectively. Hoek et al. [19] presented a proposed quantification of the GSI chart based on the rock quality designation (RQD) and the joint condition rating of the RMR system. Data from four different rock masses were used to extend the case history which was proposed by Hoek et al. [19] and provided by Bertuzzi et al., [20] and the good correlation between the GSI qualitatively assessed from the standard GSI chart and the quantified GSI was found. Han et al. [21], Poulsen et al. [22], and Wang et al. [23] also suggested different methods focusing on quantifying the GSI chart to facilitate the use of the system to determine the strength and deformation parameters of the rock mass. Morelli [24] analyzed different GSI calculation methods based on the Monte–Carlo simulations. He found the highest and the lowest GSI values from the equations which applied the conventional RMR1989 values and the RMI method, respectively. Iran is located on the Alpide earthquake belt, in the active collision zone between the Eurasian and Arabian plates. This issue makes Iran a country that suffers from geotechnical seismic hazards associated with frequent destructive earthquakes. Also, according to the rapid growth of population and demands for construction lifelines, the risk assessment studies which should be carried out in order to reduce the probable damage is necessary [25]. A principal cause of earthquake damage is landsliding, and the ability to predict earthquake-triggering landslide displacements is important for many types of seismic hazard analyses and for the design of engineered slopes [26].

Sari [27] found similar results. The digital face mapping, as a practical tool to characterize rock masses, was been investigated by [2729] which could significantly reduce the time required in the field. Hong et al. [29] proposed a method to determine the GSI quantitatively using photographic images of an in situ jointed rock mass with an image processing technology, fractal theory, and artificial neural network (ANN). Bozorgzadeh et al. [30] and Contreras et al. [31] used Bayesian statistics to quantify the uncertainty of intact rock strength. Hoek and Brown [32] introduced relatively few fundamental changes to demonstrate practical applications of the criterion and the GSI system. Day et al. [33] presented a modified GSI chart and a new Composite Geological Strength Index (CGSI) methodology to combine multiple suites of rock mass structure using a weighted harmonic average to evaluate complex rock masses. Vásárhelyi and Bögöly [34] presented a new method for calculating the GSI value using the “integral-geometric method.” It provided another GSI value calculation method that broadened the determination range of the GSI in case of poor rock mass. Hussian et al. [35] presented the review of the 19 years of research studies conducted by different researchers on the GSI. Although the GSI tool is applicable and capable of estimating the shear strength parameters of the rock mass, estimating regional scale shear strength parameters with thousands of slopes, it tends to be costly and almost impossible.

For helping the planners in selection of suitable locations to implement development projects, a landslide hazard zonation map has been produced for the Golmakan Watershed as part of Binaloud northern hillsides (northeast of Iran) [36]. In this study, we intended to develop an effective empirical regression equation to estimate GSI in Alborz mountains range in Iran, based on the generalized Hoek–Brown failure criterion and effective parameters which were obtained from geological maps, published charts and tables. In this regard, Gorgan geological quadrangle was selected as the main study area and Roudbar geological quadrangle was chosen for the verification of the model. The two regions are located at highly active seismic zone and are very susceptible to landsides, which can be triggered by any significant earthquake in the north of Iran. This model can be used to prepare initial data for the applications in the primary landslide susceptibility or hazard analysis of the entire Alborz zone at the regional scale.

2. Material and Methods

2.1. Study Area

In quantitative techniques, prediction for landslide susceptibility is based on the actual realistic data and interpretations. Further, the quantitative techniques also overcome the subjectivity of qualitative approaches [37]. From strong ground motion perspective, the Zagros thrust fault zone and Alborz and Central Iran zones are two tectonic provinces of the Iranian territory as shown in Figure 1 [26]. The study area (Gorgan geological quadrangle) with an area of 2400 km2 is located between 54° and 54° 30’ longitudes and 36° 30’ and 37° latitudes (Figure 2) at eastern Alborz tectonic zone. The area is distinguished by two major geomorphological units including the Alborz mountainous terrains in the southern parts and the Gorgan plain in northern regions. These two areas are separated by late Mesozoic Khazar (Caspian) fault. The formation of northern areas dates back to the late Tertiary as well as early Quaternary tectonics activities, causing a rapid subsidence in the north part of Khazar fault filled by a thick sequence of marine and continental deposits. These deposits were covered by aeolian deposits (known as Gorgan Loess soils) with a thickness of 5–70 meters, which were deposited during climate changes due to frozen and melting periods of Ice Age in the late Quaternary. Mountainous areas with plenty of forests, comprising the two Alborz and Gorgan-Rasht structural subzone, are separated by the North Alborz fault. The compressional tectonic regime of the region resulted in faulting and folding with overall E-W to NE-SW trends. Gorgan (Khazar) and North Alborz faults as well as Jahan-Nama syncline and Chahar-bagh anticline are some of the major structural features around the study area. A wide range of geological strata from the Palaeozoic to Quaternary periods outcrop in the study area. Metamorphic schists known as Gorgan schist with a probable age of late Precambrian is known as the oldest rock in the region [40]. Digital geologic maps of the quadrangle formed the basis of the study to assign material properties throughout the area.

2.2. GSI Estimating

As pointed out by Hoek, the GSI classification system is applicable to intact or heavily jointed rock masses. It is important that the Hoek–Brown failure criterion is widely accepted for rock masses which are assumed to be isotropic. In other words, the behavior of the rock mass would not depend on the direction of the applied loads. Therefore, the slopes in which failure surface are imposed by singular discontinuities are highly anisotropic and GSI system is not applicable [41]. When the failure plane passes through several zones, the GSI values require special judgment and the mean values may not be appropriate. A systematic study was conducted to analyze the nature and behavior of rock masses in the study area. For this purpose, firstly an 8-day field investigation into 320 geological stations was carried out by the leading authors and staff from Gorgan quadrangle. Based on the mapping of the rock exposures, all the data were collected by visual chart assessments at the scale suitable for slopes. The sampling stations were identified based on good lithological exposures and the condition of slope stability. The collected data included rock mass structure, rock type, joint condition, joint roughness, and hydrological condition. The difference image of GIS-derived landslide susceptibility zonation maps prepared for pre- and post-Chamoli earthquake shows the effect of seismic shaking on the occurrence of landslides in the Garhwal Himalaya. An attempt has been made to incorporate seismic shaking parameters in terms of peak ground acceleration with other static landslide causative factors to produce landslide susceptibility zonation map in geographic information system environment [42]. The GSI values were estimated qualitatively based on thorough geological visual field observations.

The rock mass was classified into four classes, i.e., blocky, very blocky, blocky disturbed, and disintegrated. The surface condition generally varies from smooth to poor, slightly weathered to highly weathered with soft clay infillings. Surface roughness data showed a wide range from rough to slickenside. From the hydrogeological point of view, all the visited slopes were dry. Three hundred twenty locations were considered for performing GSI calculations as shown in Figure 2. The GSI evaluation showed the typical diagonal trend from top-left to bottom-right that depicted decreasing rock mass quality. It also showed that GSI values ranged from 25 to 80, i.e., from the crushed rock to almost intact rock (Figure 1). Field photographs were taken at all slopes. The reevaluation of photographs revealed a structurally controlled behavior in some stations. Therefore, the GSI system was inappropriate. For this reason, the corresponding data of 26 stations were ignored and 294 stations were considered for representing the area. Figure 3 shows the GSI measuring for two types of rock structure.

2.3. Data Used

The main data set used in the present study was extracted from Gorgan Geological map (1 : 100,000), digital elevation model (DEM) of the studied region, and Hoek–Brown criterion table. The medium-sized geological map shows a large amount of geological information, some of which can be used to estimate the GSI. The database used in this study included lithology, distance from faults and thrusts, unit weight (γ), uniaxial compressive strength (σci), and the values of constant mί in intact rock. These data were divided into two groups. The first group consisted of the slope aspect, slope height, distance from faults, and thrust, which were extracted from the geological maps and digital elevation map (DEM). The second group included the unit weight, uniaxial compressive strength (σci), and constant mί of the intact rock which were selected from published tables and charts as discussed in the following sections.

2.4. The First Group of Data

This group of data contains the slope aspect, slope height, distance from faults, and distance from thrusts. The Gorgan 1 : 100,000 geological map (Figure 4) and Digital Elevation Model (DEM) of the study area with a resolution of 10 meters were applied for the preparation of the database. The geological map was used as a basis to extract the material properties and geological structure throughout the area. The existing rock types at this area included sandstone, shale, limestone, dolomite, marl, schist, marly limestone, conglomerate, monzodiorite, dolomitic limestone, and silty shale (Figure 4). A Digital Elevation Model (DEM) of the study area was prepared based on the digital elevation contours with intervals of 10 m. All maps were obtained from the Geological Survey and Mineral Exploration of Iran. The slope aspect is one of the basic parameters capable of influencing landslide occurrences [44]. According to the previous studies, northern slopes (N) tend to have more landslide susceptibility than southern (S) ones because of higher soil moisture and thicker soil coverage [45]. Mirsanei and Mahdavifar [12] studied 4143 landslides in Iran. They reported 1433 landslides in Alborz Mountains, the majority of which have occurred on the north and northwest aspect. In general, it seems that GSI values in northern slopes are fewer than other slope aspects due to the moisture effect on rock mass weathering and joints surface conditions.

In this study, the slope aspect was divided into eight classes: North (N), Northwest (NW), West (W), Southwest (SW), South (S), Southeast (SE), East (E), and Northeast (NE). Table 1 shows the weighting factors assigned to the different slope groups, based on the numerical scale proposed by Mirsanei et al. [12]. Outcrops were considered as a valuable source of data but they might be influenced by surface relaxation, weathering, and the alteration of rock mass components. Larger rock masses had lower overall strength compared to smaller rock masses because of scale effect. The larger rock mass involved the greater number of potential failure paths and it showed that rock masses had the ability to find failure paths with least resistance. Figure 5 shows the scale effect on rock mass compressive strength [28]. In this regard, the total-slope-height is defined as the elevation difference between the slope top and the slope toe.

2.5. The Second Group of Data-Extracted Form Charts and Tables

The second group of data comprises of the unit weight, the uniaxial compressive strength (σci), and the petrographic constant mί of the intact rock considered in Hoek–Brown failure criterion. The uniaxial compressive strength σci and the material constant mi are determined by laboratory testing or estimated from published tables. Due to the lack of experimental data, σci of the intact rock and mί were directly selected from Table 2 and Table 3, based on Hoek–Brown’s recommendation. Table 2 was the source for selecting σci as well as the updated values of mί available in Hoek et al. [46] (Table 3) or in the RocLab program (2007). The minimum value of σci and mί was assumed in order to prevent the removal of the sensitive slopes. Unit weights were selected from published tables or similar data like Barton and Choubey [47] due to the lack of experimental data.

2.6. The Regression Model for Predicting GIS

Determining the rock mass shear strength parameters including internal friction angle (φ) and cohesion (c) according to Mohr–Coulomb criterion is an important step in performing landslide susceptibility analysis or hazard mapping based on Newmark’s simplified displacement method. The method to determine the strength parameters depends on the size of the study area. Hoek et al. [48] presented equations to find equivalent Mohr–Coulomb parameters predicted by the Hoek–Brown failure criterion. Hence, by estimating the GSI values, the equivalent c and φ could be obtained. Due to the limitation of the common methods such as the laboratory and field testing and failure criteria for the regional scale study, a regression model was applied to develop an empirical formula to predict the GSI. In the first step, a linear regression model was considered to determine the model coefficients. As mentioned earlier, the database contained seven input parameters. In order to investigate the effect of predictors on GSI and also to propose a comprehensive model for determining GSI, the regression (MR) analysis was performed. The aim of this method was to develop a linear regression equation for approximating science and engineering problems. When there is more than one input model, MR can be performed to obtain the best-fit equation. By using a constructed database consisting of 294 geological stations datasets in which the distance to fault (FD), uniaxial strength of intact rock (σci), rock constant (mί), slope height (H), slope aspect, and distance to fold axis (DFA) were considered as model predictors to estimate GSI, IBM SPSS Statistics 20 program was employed and the MR model was developed. Table 4 shows the summary of the model.

3. Results and Evaluation

The value of R2 is 0.86 and the adjusted R2 is 0.73 showing that seven predictors entered in the regression analysis account for 73% of the variations of GSI (Table 4) and this will fit the data at a very high level of statistical significance. Table 5 shows that the F-value is equal to 89.4 significant to a level of 5%, which is much greater than 1. Therefore, this indicates a linear relationship between the variables. The following parameters were found to be significant for GSI: distance to fault, distance to fold axis, and σci. The parameters that were found insignificant are mi, slope height, slope aspect, and unit weight. Table 5 shows the computation for the model. The estimation regression model for the GSI is presented as follows:where FD is the distance to fault, σci is the uniaxial strength of intact rock, and the last term is DFA, which is the distance to the fold axis of the model.

3.1. Testing the Assumptions of the Regression Analysis

The linear regression model is based on four assumptions. These postulate the properties that the variables should have in the population. The regression model only provides proper inference if the assumptions are held true (although the model was robust to mild violations of these assumptions). The reliability and validity of the model were considered by checking the model assumptions of the null hypothesis, irrespective of the residuals, linearity, and normality.

3.1.1. F-Test

The F-test is used in regression analysis to test the hypothesis that all model parameters are zero versus the alternative that at least one slope coefficient is nonzero. The Hypothesis test determines whether there is a relationship between the response and the predictor. When there is no relationship between the response and any of the predictors, the model will not explain much of the variation in the response. The Mean Square Model and Mean Square Error will be approximately the same, and the F Ratio will be close to 1. On the other hand, if the alternative hypothesis is true, at least one coefficient is nonzero. The Mean Square Model will be greater than the Mean Square Error, and the F Ratio will be greater than 1. According to Table 6, the F Ratio was 89.437 and it was concluded that at least one term in our model was significant.

3.1.2. The Durbin–Watson Test

The Durbin–Watson (DW) statistic is a test for autocorrelation in the residuals from a statistical regression analysis. The Durbin–Watson statistic will always have a value between 0 and 4. A value of 2.0 means that there is no autocorrelation detected in the sample. Values from 0 to less than 2 indicate positive autocorrelation and values from 2 to 4 indicate negative autocorrelation. As shown in Table 4, d = 1.624 is the critical value in the range of 1.5 < d < 2.5. Therefore, it could be assumed that there was no first-order linear autocorrelation in the multiple linear regression data. In addition, the errors were independent, and as a result, the utility of the regression model was confirmed.

3.1.3. The Multicollinearity

Multicollinearity is known to be undesirable when one independent variable is a linear function of other independent variables. According to the last column of Table 5, the value of VIF in all independent variables was less than 5 (VIF<5). Accordingly, there was no multicollinearity between independent variables; hence it validated the fitted model.

3.1.4. The Normality

Normality is defined as the normal distribution of residuals in predicted responses with an average of zero. Figure 6 illustrates normality checks in this study. For the standardized residual histogram to appear normal, a fitted normal distribution aid was considered. The average value, presented at the right side of the diagram, was very small (close to zero) and the standard deviation was approximately equal to unity. Figure 7 shows the P-P plot, which plots variable cumulative proportions against the cumulative proportions of any number of test distributions. Probability plots are generally used to determine whether the distribution of a variable matches a given distribution. If the selected variable matches the test distribution, the points would cluster around a straight line. The estimation error is also calculated for the comparison. Figure 8 shows the error values based on

and are the measured and predicted values of GSI, respectively. The estimation error is between -15 and 15, while the average error is between -5 and 5.

3.2. The Verification of Formula and Results

In this section, the applicability of the regression equation in Alborz zone using the Roudbar geological quadrangle data set was verified. This data set was prepared based on (1) and includes σci, the distance to fault, and the distance to the fold axis. The testing area (Roudbar geological quadrangle) is 2400 km2, located at 49° and 49° 30’ longitudes and 36° 15’ and 37° latitudes (Figure 9).

It faces Alborz zone and the Caspian Sea with a humid climate which might deeply affect the occurrence of landslides, and it is also mountainous (the height differs from 150 to 2,800 meters above the sea level). All of these properties along with the seismic potential of the region make it a landslide-prone area. The study area tends to comprise a wide range of geological materials, from the metamorphic rocks of Precambrian to Quartz sediments. In Alborz, the red sandstone deposits can be seen almost everywhere, from the Infracambrian, with no changes in sedimentation and geological characteristics, to the Palaeozoic. There are no Silurian and Carboniferous, representing two top volume bulges; i.e., Triassic rocks are left uncomfortably overlying Permian deposits, with the other being observed at the base of Shemshak formation (Jurassic) [49]. 1 : 100,000 geological map and Digital Elevation Model (DEM) of the study area with the resolution of 10 meters were applied (Figure 10). These maps were obtained from the Central Geological Survey and mineral exploration of Iran.

The field data were collected from 300 slopes of varied geological and slope stability conditions. The selected slopes presented a variety of rock types having various discontinuities. Surface condition generally showed a wide range of quality with different infilling cases. Three-hundred locations were considered for performing GSI calculations as shown in Figure 11. Equation (1) was applied for this area and all locations mapped GSI were compared with predicted values. Figure 12 shows that the GSI difference values are between -30 and +30 and the mean error is between −5 and +5. This result shows that (1) can estimate GSI values with an average error percentage of about 3%.

3.3. Discussion

The ability to predict landslide displacements is important for many types of seismic hazard or susceptibility analysis and Newmark’s approach provides a useful means to predict landslide displacement. This method requires knowing the static factor of safety and the landslide geometry. Determining the representative shear strength values for geological formations is necessary to calculate the factor of safety. Although the GSI system could be considered as a more appropriate and rapid method in comparison with other methods, this method needs huge field studies. In this study, an empirical regression equation was developed among some geological and geomorphological parameters, extricable from geological maps and DEM of an area and GSI values of rock masses that can reduce large field works. The parameters used to generate the model were obtained from digital geological map, DEM, and published charts and tables. For the sample area, Gorgan geological quadrangle, the database, which was prepared, contained 7 parameters. The developed model revealed that only 3 of 7 parameters made significant contribution to the model, while 4 of them made insignificant contribution; therefore, they were removed from the model. In order to test the model, another area (Roudbar geological quadrangle) was considered. Figure 11 shows the comparison between the results of GSI obtained from empirical equation and the field work. All the results indicated the GSI values predicted by the empirical equation were generally in good agreement with field data. On the other hand, the dispersion of the predicted GSI does not follow any particular pattern. Therefore, it is not clear which predictor variable has a significant effect on the dispersion of the GSI. There are many simplified assumptions for regional assessment of input data, in spite of the mentioned limitations. It should be noted that this method provides a rapid estimation of shear strength parameters and can be used for emergency responses.

4. Conclusion

Landslides are one of the most damaging geotechnical hazards associated with earthquakes. In addition, the prediction of triggered areas after an eventual earthquake is an important issue for engineers as well as risk managers. The shear strength parameters of geological formation are important elements in such types of analyses. In the present study, a simple method based on the geological map and published tables were used to develop a regression equation between GSI value and easily accessible parameters for a study area. For this study, an area was selected as a pilot area (the Gorgan geological quadrangle in the northeast of Alborz region). The regression relationship was extracted based on observation in 294 natural or excavated slopes, their GSI identification, and the linear correlation determination between these GSI values and some geological and geomorphological parameters. The regression model presented in (1) is well constrained (R2 = 74%) and predicts GSI in terms of (1) distance from the fault (FD), (2) the uniaxial strength of the intact rock (σci, determined from published table for different rock types), and (3) the distance from fold axis (DFA). The approach used is so simple that results can easily be updated. In the next steps, another area (the Rudbar geological quadrangle in the west of Alborz) was selected as the validation area and the evaluation of predictive ability of regression relationship. This validation step showed that the accuracy of this model was in an acceptable range. The average error of this model ranges between ±5. To examine the supplication of the model in other regions, Roudbar region at the south of Gilan province with geological and seismic similarity was considered to predict the GSI and compare it with the observed data [23, 50, 51].

Data Availability

The data used to support the findings of the study can be obtained from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.