Using multivariate adaptive regression splines to develop relationship between rock quality designation and permeability

Mohsin Usman Qureshi, Zafar Mahmood, Ali Murtaza Rasool

a Faculty of Engineering, Sohar University, Sohar, Oman

b University of Buraimi, Al Buraimi, Oman

c National Engineering Services Pakistan, Lahore, 54700, Pakistan

Keywords:In situ permeability Limestone Sandstone Lugeon Rock quality designation (RQD)Multivariate adaptive regression splines(MARS)

ABSTRACT The assessment of in situ permeability of rock mass is challenging for large-scale projects such as reservoirs created by dams,where water tightness issues are of prime importance.The in situ permeability is strongly related to the frequency and distribution of discontinuities in the rock mass and quantified by rock quality designation (RQD). This paper analyzes the data of hydraulic conductivity and discontinuities sampled at different depths during the borehole investigations in the limestone and sandstone formations for the construction of hydraulic structures in Oman.Cores recovered from boreholes provide RQD data, and in situ Lugeon tests elucidate the permeability. A modern technique of multivariate adaptive regression splines(MARS)assisted in correlating permeability and RQD along with the depth.In situ permeability shows a declining trend with increasing RQD, and the depth of investigation is within 50 m.This type of relationship can be developed based on detailed initial investigations at the site where the hydraulic conductivity of discontinuous rocks is required to be delineated. The relationship can approximate the permeability by only measuring the RQD in later investigations on the same site, thus saving the time and cost of the site investigations.The applicability of the relationship developed in this study to another location requires a lithological similarity of the rock mass that can be verified through preliminary investigation at the site.

1. Introduction

The success of geotechnical projects having water tightness issues is directly related to the accurate assessment of in situ permeability with cost-effectiveness. The geotechnical practice employs the widely accepted and most reliable Lugeon water pressure tests to delineate in situ permeability in rock mass(Houlsby, 1976). Currently in practice, the techniques used to delineate the hydraulic properties are expensive and timeconsuming (Long and Witherspoon, 1985). Since the flow occurs through the discontinuous rock mass, their properties should be used to characterize the hydraulic conductivity (Witherspoon and Gale, 1982). Zhang (2013) highlighted the effects of discontinuity frequency in fractured rock on the anisotropy of permeability.

In the past, many researchers (Magnusson and Duran, 1984;Nappi et al., 2005; Hamm et al., 2007; Jiang et al., 2009; Qureshi et al., 2013, 2014; Angulo et al., 2011; Shahbazi et al., 2020) reported the estimation of hydraulic conductivity of rock mass employing discontinuity data. A comparative study (Magnusson and Duran, 1984) between the electrical resistivity and hydraulic properties of fractured core logs concluded that the crystalline rock’s micro-fractures mainly transport the electric current and groundwater.The study also concluded that"a few discrete fractures govern the hydraulic conductivity and the variation between two methods is due to the effect of restricted water flow on few channels with small apertures, while electric surface conduction was possible through fractures with tiny apertures"(Magnusson and Duran,1984).In another study, Hamm et al. (2007) used an acoustic televiewer and recovered granite cores to study the relationship between the hydraulic conductivity and fracture properties estimated from packer tests and borehole data in fractured granite,and concluded that the hydraulic conductivity is more strongly related to the aperture of fracture compared to its frequency.Angulo et al.(2011)investigated the borehole data, including electrical resistivity,fracture frequency, and low-pressure permeability test (LPT) results. The investigation depths were 100 m and 120 m in two boreholes drilled in karst formation, mainly composed of limestone. Their study concluded that hydraulic conductivity shows a more significant relationship with electrical resistivity than discontinuities. Nappi et al. (2005) proposed a relationship between the hydraulic conductivity and depth of rock mass for the initial characterization of fractured media’s hydraulic properties.Outcrop measurements and Lugeon tests were employed to investigate the hydraulic characteristics of sandstone.

In the case of a site with no history of tectonic events and uniform lithology, a depth-wise decrease in the discontinuity’s frequency could be related to decreasing mean hydraulic conductivity and an increase in the mean rock quality designation (RQD) (Jiang et al.,2009). According to Zhang (2013), the hydraulic conductivity of discontinuous and intact rock masses is affected by geostatic stress.Usually, permeability decreases with depth, with a high rate of change near the surface,but the rate decreases with depth(Shahbazi et al.,2020).Various numerical modeling approaches concluded that the permeability is directly related to the shear stress (Xiong et al.,2011)and inversely related to normal stress(Chen et al.,2015).

Brown (1981) proposed RQD, a standardized method to quantitatively describe the discontinuity characteristics of rock mass.RQD is easy to be obtained from the visual observation of core recovered from a borehole. In the existing set of equations, the permeability of discontinuous rocks is correlated with rock mass rating (RMR) (El-Naqa, 2001), RQD (El-Naqa, 2001; Qureshi et al.,2014), and geological strength index (GSI) (?ge, 2017) which shows a strong impact of discontinuity frequency on hydraulic conductivity. RQD is the most widely accepted among the considered rock mass indices and can be easily obtained from the recovered cores during borehole drilling. A linear regression relationship between the mean RQD and depth was reported by Jiang et al. (2009), and more recently, empirical equations to assess the hydraulic conductivity of rock mass with respect to depth were presented by Chen et al.(2018)and Piscopo et al.(2018).Although the studies showed the importance of the depth factor,the authors also believed that the effects of geostatic stresses are already included in the in situ hydraulic properties of the discontinuous rock mass. The empirical correlations reported by El-Naqa (2001)and Jiang et al. (2009) are shown in Eqs. (1) and (2), respectively.

where KLuis the permeability measured in Lugeon(Lu),and K is the permeability measured in m/d.

The present research is an extension of Qureshi et al.(2014).The permeability is correlated to the RQD of discontinuous sedimentary rocks by utilizing the hydraulic and fracture properties obtained during the site investigations for Wadi Dayqah dam, Oman (Prisk et al., 2009). The initial relationship developed by Qureshi et al.(2013), applicable to the local sedimentary rocks having RQD >25%, is given as

where K is measured in cm/s.

Later, a modified relationship was reported by Qureshi et al.(2014) as

In the current study, the multivariate adaptive regression splines (MARS) technique, developed by Friedman (1991), was utilized to show the relationship between the permeability and RQD and depth. MARS can produce simple-to-interpret models by extracting complex data mapping in high-dimensional data.Another advantage of this method is to express the relationship mathematically.Furthermore,the MARS can address the analysis of geotechnical systems (Zhang and Goh, 2013; Zhang et al., 2016).Similar modern techniques of analysis have been employed in the field of geotechnical engineering, i.e. evaluation of liquefaction potential (Zhang et al., 2015, 2016), stability assessment of excavations(Goh et al.,2017),inverse analysis of soil and wall properties in braced excavation (Zhang et al., 2017), determination of tunnel related earth pressure balance(Goh et al.,2018),reliability analysis of pile foundations (Kumar and Samui, 2019), assessment of pile drivability (Zhang et al., 2021a), determination of horizontal wall defection envelope for braced excavations in clay (Zhang et al.,2019), prediction of tunnel displacement due to excavation using MARS (Zheng et al., 2020), prediction of undrained shear strength(Zhang et al., 2021b), and probabilistic stability analysis of earth dam slope (Wang et al., 2020). Samui (2020) reviewed the application of different artificial intelligence techniques in geotechnical engineering, including MARS. Considering the versatility of the MARS algorithm for application in geotechnical engineering,in the current study,the relationship developed by Qureshi et al.(2014)is further enhanced by employing the MARS algorithm and including the factor of depth.

2. Methodology

The in situ data of RQD and permeability recorded within a depth of 50 m are utilized for developing a significant relationship using MARS. A brief description of the methods used and the material properties are described in the following sections.

2.1. Multivariate adaptive regression splines (MARS)

MARS method is based on a statistical approach for fitting the data of dependent and independent variables by a series of splines of different gradients (Friedman, 1991). The endpoints of the splines(knots)mark the end of one set of data and the beginning of another, resulting in piecewise curves called basis functions or hinge functions. In the present study, analysis of the MARS model was performed using earth package (Milborrow, 2018) in R Studio Version 4.0.0 (R Core Team,2020).

Let y be the response variable and X = (X1,…, XM)be a matrix of M input variables. The continuous response model is shown as follows:

where e is the error. The function f is approximated by applying basis functions. Basis functions are splines, including piecewise linear functions of the form (x-α)+, where the knot occurs at x = α. The term (x-α)+means that only a positive part is used;otherwise, it is zero (Fig.1).

Formally basis function is expressed as

Fig.1. One-dimensional linear spline with knot (or hinge) at α = 0.5.

The MARS model is expressed as (Zhang et al., 2016):

where β0is a constant,βiis the coefficient of the ith hinge function,hi(X)is a hinge function,and n is the number of hinge functions of the model.

Model subsets are compared using generalized cross-validation(GCV). The GCV is a form of regularization that trades off the goodness-of-fit against the model complexity.The GCV of a model is calculated as follows (Hastie et al., 2009):

where M is the number of basis functions, and d is the penalizing parameter.The optimal value of d usually falls in the range of 2 ≤d≤4, and generally d = 3 is used (Friedman,1991).

The residuals are the difference between the values f(x) predicted by the model and corresponding response values y. The residual sum of squares (RSS) is the sum of the squared values of residuals:

The coefficient of determination R2or RSq measures the variance in response variable y that can be predicted using f(x). The RSq standardizes the residual sum of squares (RSS) and is defined as follows:

Generalized R2or GRSq is the generalization performance of the model estimated using the MARS algorithm. GRSq is defined as follows (Milborrow, 2021):

GCV is important for MARS because it is used to evaluate model subsets in the backward pass.

2.2. Rock quality designation (RQD)

Deere (1964) first introduced the RQD index, a well-known parameter, to characterize the distribution of discontinuity in the rock mass.RQD is reported as the percentage of the rock core length(l)comprising intact lengths of at least 10 cm(li).The mathematical expression for calculating the RQD is given as

Later,Deere and Deere(2009)signified the versatility of RQD in the practice of hydraulic structures,and it was standardized by the International Society for Rock Mechanics and Rock Engineering(Brown, 1981). Visual observation of the rock cores recovered during the drilling of boreholes quantifies the RQD.

2.3. Lugeon water pressure test

The Lugeon tests delineate the permeability in the borehole at the specified depth intervals during the drilling operation, as shown in Fig.2.The Lugeon water pressure test(Houlsby,1990)is a constant head type test in which water at constant pressure is injected into the borehole through a perforated pipe bound by a double packer system (Fig. 2), and the discharge is measured. The maximum test pressure, Pmax, which should not exceed geostatic stress expected at the center of the respective section,is calculated as 1 psi (6.897 kPa) per foot (0.3048 m) depth of the borehole.During the test,water is injected at five different pressures,i.e.50%,75%,100%,75%,and 50%of Pmax,each being maintained for 10 min(Houlsby,1990), as shown in Fig. 3.

Fig.2. Schematic diagram of field tests-Lugeon test and RQD(Deere,1964;Houlsby,1990).

Fig. 3. Different water pressures during Lugeon test (Houlsby,1990).

The Lugeon value (Lu) at each pressure is calculated using the following equation:

where α equals 1 in the international system of units, q is the discharge,P0is the reference pressure equivalent to 1 MPa,L is the section length (Fig. 2), and P is the applied pressure at any stage.The procedure given by Houlsby (1976) is employed for the selection of Lugeon value which is interpreted from the flow pattern during the test. The Lugeon value was correlated by Fell et al.(2005), and according to their study,1 Lu = 1.3 ×10-5cm/s.

3. Results and discussion

The current study analyzed the results of 61 Lugeon water pressure tests performed within a depth of 50 m in the discontinuous sedimentary rock formations of the Wadi Dayqah dam project site (Prisk et al., 2009). The project site is located in a topographically steep area underlain by a sequence of sedimentary and metasedimentary rocks. It lies on a faulted tertiary block, and rocks are fractured and jointed. Despite their highly crushed nature, the permeability of the rocks is relatively small, and the rocks are in a compacted state due to the thrusting of faults.Within the depth of 50 m, the lithology of limestone was described as grayish to brownish yellow-bioclastic, nodular, and micritic, and the sandstone as whitish-gray,poorly cemented,fine-grained,massive,and friable. The limestone has been affected by solution and is characterized by the presence of solution enlarged joints (Qureshi et al.,2014). The current study analyzed a total number of 61 Lugeon tests, 37 were performed in limestone and 24 in sandstone formations. The scatter of data is shown in Fig. 4. The mean depthwise curve is also plotted for RQD and permeability, which shows a general trend of increase in RQD and decrease in permeability with depth. The permeability-depth equation proposed by Chen et al. (2018) and the RDQ-depth equation obtained by Jiang et al.

(2009) are also plotted in Fig. 4a and b, respectively. The increase in RQD and decrease in permeability with depth indicate the possibility of an interrelation between RQD and hydraulic conductivity.

Fig. 4. The scatter plots of (a) permeability (K) and (b) RQD versus depth.

Fig. 5. Frequency distribution of (a) permeability (K) and (b) RQD obtained from 61 Lugeon tests.

The frequency distribution of permeability and RQD recorded in the field for all tests is shown in Fig. 5. The sandstone exhibits higher permeability and lower RQD as compared to limestone.The majority of permeability data recorded range between 10-5cm/s and 10-3cm/s for sandstone and between 10-6cm/s and 10-3cm/s for limestone.

K-fold cross-validation is a statistical method to estimate the performance of machine learning models.Typical choice of K-folds is 5 or 10(Kohavi,1995).While building a machine learning model,cross-validation helps to ensure that the model fits the data accurately without overfitting. The noisy estimate of model performance may be the result of a single run of K-fold cross-validation.The variation in results may occur due to different splits of data.Repeated K-fold cross-validation improves the performance of machine learning model, as cross-validation process is repeated multiple times. In this study, 3 times repeated 10-fold crossvalidation is used to obtain the model. The model is trained on 9 folds of the data and tested on the 10th fold, and the process is repeated 10 times by changing the test fold. The process of crossvalidation is repeated 3 times. The hyperparameters of the model are shown in Table 1.

Selection of optimum number of terms of MARS model is shown in Fig.6.Thin gray lines show the out-of-fold RSq for each fold.The thick red line is the mean of the gray lines.The thick black line is the GRSq line for the full model, ignoring cross-validation. The dotted blue line is the RSq line for the full model,ignoring cross-validation.RSq and GRSq are plotted against the number of terms used in the model selection. The RSq keeps increasing with increasing terms(i.e. with increasing model complexity), however, GRSq starts decreasing with model complexity.The final model with 5 terms is selected when the mean out-of-fold RSq is maximum.

Table 1 The hyperparameters for MARS model.

Fig. 6. Model selection in a 3 times repeated 10-fold cross-validation.

Fig.7 is the residual versus fitted graph that shows the residual for each of the fitted response values. Residuals show constant variance, i.e. the residuals remain evenly spread-out or homoscedastic as the fitted values increase. However, constant variance is not as important in the MARS model as in the linear models(Milborrow, 2021). The red curve is a lowess (locally weighted scatterplot smoothing) fit. It can be recalled as fancy moving averages.Lowess fit is close to zero and does not show any divergence at low or high fitted values.

Fig. 8 shows the cumulative distribution of absolute values of the residuals.The gray-colored curve is the distribution of absolute residuals, while the black-colored step-shaped curve is the cumulative distribution of absolute residuals. The bell-shaped distribution of absolute residuals translates to an approximate S-shaped cumulative distribution curve (which is not applicable here as the distribution of absolute residuals is not bell-shaped). The median absolute residual is 0.746(vertical gray line for 50%).The 95%of the absolute values of residuals is less than 2.17. Therefore, in the training model, at 95% of the time, the predicted value is within ln2.17 or 0.77 times of the observed value.

Fig. 7. Residual versus fitted graph of the final model.

Fig. 8. Cumulative distribution of absolute values of residuals of the final model.

Fig. 9. Quantile-quantile plot of the final model.

Fig.9 shows a quantile-quantile plot in which the distribution of residuals to normal distribution is compared. Since normal and residual quantities lie on the black dotted line, the residuals are distributed normally.This plot is useful in figuring out the outlying residuals or outliers.The distribution of residuals(marked“actual”in Fig. 9) matches with the normal distribution (marked “normal”in Fig. 9). Since all data lie on the black dotted line, there is no divergence of residuals (also evident in Fig. 7).

The model’s performance is determined by RSq(or R2)(Eq.(12)),root mean squared error(RMSE),mean absolute error(MAE),mean error (ME), and relative absolute error (RAE) according to the following equations:

where f(x)is the model predicted value.The RSq(Eq.(12))and GRSq(Eq.(13))of the final model are 0.766 and 0.6631,respectively.The performance metric of the final model is shown in Table 2.

The final MARS model equation is shown as follows:

where K is in cm/s,and D is the depth in m.The model developed in this study shown in Eq. (20), the relationship derived by Qureshi et al. (2014) and shown as Eq. (4), and the examined data for this study are plotted in Fig. 10. The MARS model best fits the data analyzed in the study and demonstrates a declining trend of permeability with the increase of RQD. The MARS model also includes the depth factor which has a prime importance for the evaluation of permeability as concluded by Chen et al. (2018) and Jiang et al. (2009).

The relationship developed using the MARS model and shown as Eq.(20)is a result of analysis of data obtained from the boreholeinvestigations and is proposed to access the permeability by sampling the RQD at the required depth for the preliminary site investigations.It is important to note that the relationship presented in Eq. (20) is based on the data obtained from 61 Lugeon tests performed within a depth of 50 m. However, it is a useful tool to make an initial assessment on the hydraulic characteristics of discontinuous rocks of the site in question. Further, investigations in the same geological setting can include the measurement of RQD only or increase the depth interval of the permeability tests. This approach can save the time and cost of conducting expensive and time-consuming in situ permeability tests reporting similar results.It can also be used for comparison purposes to achieve confidence in the site having similar geological properties of rock as reported in this study. A good agreement between the permeability values obtained by MARS model (Eq. (20)) and actual recorded data is shown in Fig.11.

Table 2Performance metrics of MARS model.

Fig.10. Plot between the permeability (K) and RQD.

Fig. 11. Comparison of the permeability (K) calculated by MARS model and actual values.

4. Conclusions

In the present study, rock discontinuity frequency is correlated with hydraulic conductivity by the inclusion of depth factor.Lugeon tests and visual observation of rock cores elucidated hydraulic conductivity and RQD, respectively, which were obtained from borehole investigation in sedimentary rocks within a depth of 50 m.Finally, the interrelation between permeability and RQD is developed, and a comparison is made between actual and estimated values from the MARS model.The main conclusions for the current study are stated below:

(1) At the investigation site, an increasing trend of RQD with depth and a decreasing trend of permeability are recognized,which are similar to the observations reported earlier.

(2) A modern data analysis technique, MARS, is employed to develop a significant relationship between permeability,RQD, and depth. The relationship can approximate the permeability (in cm/s) by inputting the recorded RQD (in %)at a certain depth (in m).

(3) The relationship developed in the study can save the cost and time for the detailed planning and execution of geotechnical investigations on projects of similar nature.Developing such a relationship using the initial data can assist in approximating permeability values by only sampling the RQD index at the required depth in the later investigations. The relationship can also assist in deciding on reducing the number of permeability tests in the same borehole by increasing the depth interval of the test if similar RQD values are recorded.

(4) The development of relationship using MARS in the current study is a useful tool for assessing the in-situ permeability of discontinuous rock formation for a particular site. Later, the relationship can assist in the detailed planning of tests and assess hydraulic properties of rock for further investigations.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors are grateful to the engineers and experts of National Engineering Services Pakistan for their guidance and support. The authors are also indebted to the Sohar University and the University of Buraimi, Oman, to support this study.