Pinnaduwa Kulatilake
 Visiting Scholar, Materials Science and Engineering
Contact
 (520) 6216064
 MINES, Rm. 141
 TUCSON, AZ 857210012
 kulatila2019@yahoo.com
Bio
No activities entered.
Interests
No activities entered.
Courses
201819 Courses

Prob+Stat Geologic Media
MNE 502 (Spring 2019)
201718 Courses

Dissertation
GEN 920 (Fall 2017) 
Rock Slope Design
MNE 529 (Fall 2017)
201617 Courses

Dissertation
GEN 920 (Spring 2017) 
Prob+Stat Geologic Media
GEN 402 (Spring 2017) 
Prob+Stat Geologic Media
MNE 402 (Spring 2017) 
Prob+Stat Geologic Media
MNE 502 (Spring 2017) 
Dissertation
GEN 920 (Fall 2016) 
Rock Slope Design
MNE 529 (Fall 2016)
201516 Courses

Blck Thry Apps Rock Mass
CE 543 (Spring 2016) 
Dissertation
GEN 920 (Spring 2016) 
Independent Study
GEN 599 (Spring 2016) 
Research
GEN 900 (Spring 2016) 
Thesis
GEN 910 (Spring 2016)
Scholarly Contributions
Journals/Publications
 Yunfeng, G. e., Kulatilake, P. H., Tang, H., & Xiong, C. (2014). Investigation of natural rock joint roughness. Computers and Geotechnics, 55, 290305.More infoAbstract: The paper provides a comprehensive review on rock joint roughness measurement and quantification procedures. Superiority of fractal based methods over JRC, statistical parameters and statistical functions in quantifying roughness is discussed in the paper. Two of the best fractal based methodologies available in the literature, the modified 2D divider and variogram methods, are used to quantify natural rock joint roughness in 3D and 2D, respectively. The capability of these two methods in accurate quantification of natural rock joint roughness is shown in the paper by applying the procedures to four natural rock joints. A good comparison has been obtained from the values obtained through the two methods. Both these methodologies have two parameters to capture the stationary roughness. The fractal dimension captures the spatial auto correlation of roughness; the other parameter captures the amplitude of roughness. Anisotropic roughness has been studied by applying two other methodologies: (a) a triangular plate methodology and (b) a light source methodology to the same four natural rock joints. A reasonably good comparison has been obtained through the results of these two methodologies. All four roughness quantification methodologies can be applied to any size of sample covering from laboratory to field scales. The results of the triangular plate and light source methodologies provided possible sliding direction values (under the gravitational loading) close to that reported in the literature for the rough discontinuity planes used in the study. © 2013 Elsevier Ltd.
 Kulatilake, P. H., & Wu, Q. (2013). Development of an orthotropic constitutive model for a jointed rock mass. 47th US Rock Mechanics / Geomechanics Symposium 2013, 1, 555562.More infoAbstract: Fracture data available for a limestone rock mass were used to build and validate a stochastic 3D fracture network model. A procedure is proposed to investigate the size effect and the REV and equivalent continuum behaviors of fracture and mechanical properties in 3D of the jointed rock mass having finite size fractures capturing the anisotropic effects arising from the fracture system. An REV size of about 25 m, which is about 8 times the mean fracture size of joint sets, was found to represent the mechanical properties of the rock mass. A number of relations are developed between the rock mass mechanical parameters and fracture tensor components in 3D. Based on the mechanical parameter values obtained in every 45 degree direction in 3D, the principal values, principal directions and tensors are developed for rock mass mechanical parameters to represent the REV size properties. An incrementally linear elastic, orthotropic constitutive model is suggested to represent the equivalent continuum prefailure mechanical behavior of the rock mass by incorporating the effect of fracture geometry by the fracture tensor components. Copyright 2013 ARMA, American Rock Mechanics Association.
 Kulatilake, P. H., & Wu, Q. (2013). Tunnel stress analyses in 3D using equivalent continuum and discontinuum procedures. 47th US Rock Mechanics / Geomechanics Symposium 2013, 4, 27742781.More infoAbstract: Deformation and stability around a tunnel in a limestone rock mass is investigated using both (a) an equivalent continuum/discontinuum model and (b) a fully discontinuum model through threedimensional stress analyses. In model (a), the Representative Elementary Volume (REV) mechanical properties estimated in a previous paper are used to represent the combined behavior of intact rock and 4 of the 5 fracture sets that had significantly smaller sizes compared to the fifth fracture set (bedding planes). To this system, the bedding plane fractures were added according to its statistical geometric properties. In model (b), all 5 fracture sets were generated according to the obtained fracture geometry statistical models. Using model (a), effect of several factors, such as REV property values, bedding plane modeling method, the average horizontal stress to the vertical stress and tunnel depth, on the geotechnical behavior around the tunnel was studied in detail. A comparison is made between the stress analyses results obtained through the aforementioned two models. Copyright 2013 ARMA, American Rock Mechanics Association.
 Kulatilake, P. H., Qiong, W. u., Zhengxing, Y. u., & Jiang, F. (2013). Investigation of stability of a tunnel in a deep coal mine in China. International Journal of Mining Science and Technology, 23(4), 579589.More infoAbstract: Stability level of tunnels that exist in an underground mine has a great influence on the safety, production and economic performance of mines. Ensuring of stability for softrock tunnels is an important task for deep coal mines located in high in situ stress conditions. Using the available information on stratigraphy, geological structures, in situ stress measurements and geomechanical properties of intact rock and discontinuity interfaces, a threedimensional numerical model was built by using 3DEC software to simulate the stress conditions around a tunnel located under high in situ stress conditions in a coal rock mass in China. Analyses were conducted for several tunnel shapes and rock support patterns. Results obtained for the distribution of failure zones, and stress and displacement fields around the tunnel were compared to select the best tunnel shape and support pattern to achieve the optimum stability conditions. © 2013 Published by Elsevier B.V.
 Nengxiong, X. u., Kulatilake, P. H., Tian, H., Xiong, W. u., Nan, Y., & Wei, T. (2013). Surface subsidence prediction for the WUTONG mine using a 3D finite difference method. Computers and Geotechnics, 48, 134145.More infoAbstract: WUTONG coal mine is adjacent to an auxiliary dam of the Yuecheng Reservoir. In this paper, mininginduced surface subsidence prediction is conducted by means of the finite difference method (FDM) to judge whether the extraction of the coal seam will have a negative impact on the dam. First, the initial values of the rock mass mechanical parameters are estimated using the available literature that relates intact rock and discontinuity properties to rock mass parameters. Then, based on available surface subsidence monitoring data on WUTONG's mined areas, the main mechanical parameters of coal and rock masses are determined by a back analysis procedure that combines an experimental design technique with numerical simulations. Finally, the surface subsidence results in the mining area are numerically predicted for four different mining scenarios (S1 through S4). Scenario S3 emerged as the best choice of these four scenarios. The predictions are: (1) the maximum surface subsidence within the mining area is 2.14. m, with the maximum settlement point located in the midwest area of the coalfield, and (2) the nearest distance from the boundary of the surface movement area to the edge of the dam foundation is 35. m. Therefore, mining the coal seam will not cause damage to the dam. © 2012 Elsevier Ltd.
 Ren, Q., Tang, H., & Kulatilake, P. H. (2013). The deformation behavior of the soft rock applied in the construction of the underground pipe work engineering. ICPTT 2013: Trenchless Technology  The Best Choice for Underground Pipeline Construction and Renewal, Proceedings of the International Conference on Pipelines and Trenchless Technology, 812817.More infoAbstract: The deformation behavior of the rock and soil is a key problem in the construction of the underground pipe work, and will directly influence the time limit and difficulty of a project. Actually for the underground rock mass, the construction of pipe work is a stress release process which makes the deformation characteristic of the rock mass totally different with other conventional constructions, such as building construction, road and bridge construction. The reason is the stress paths of the two kinds of engineering are different; one is download stress path and the other is upload stress path, so the mechanical property and deformation characteristic of a kind of soft rock under both download and upload stress paths in the Three Gorges area were studied by using the elastoplastic modeling method. Through the experiments of lamellar marlstone of Badong formation, the manner of the interaction between plastic volume stain and shear strain, the direct effect of the plastic volume stain on shear resistance were investigated, and the conclusion that stress path is exactly the comprehensive performance form of the interaction between plastic volume stain and shear strain which has been proved by the trixal experiments of soft rock. Finally, the deformation behavior of marlstone was used to calculate the settlement in the numerical simulation of the underground pipe work's construction. © ASCE 2013.
 H., P. (2012). Preface: Special issue on selected topics in rock mechanics and rock engineering. Geotechnical and Geological Engineering, 30(3), 523524.
 Kulatilake, P. H. (2012). Preface: Special Issue on Selected Topics in Rock Mechanics and Rock Engineering. Geotechnical and Geological Engineering, 12.
 Kulatilake, P. H., Hudaverdi, T., & Qiong, W. u. (2012). New prediction models for mean particle size in rock blast fragmentation. Geotechnical and Geological Engineering, 30(3), 665684.More infoAbstract: The paper refers the reader to a blast data base developed in a previous study. The data base consists of blast design parameters, explosive parameters, modulus of elasticity and in situ block size. A hierarchical cluster analysis was used to separate the blast data into two different groups of similarity based on the intact rock stiffness. The group memberships were confirmed by the discriminant analysis. A part of this blast data was used to train a singlehidden layer back propagation neural network model to predict mean particle size resulting from blast fragmentation for each of the obtained similarity groups. The mean particle size was considered to be a function of seven independent parameters. An extensive analysis was performed to estimate the optimum value for the number of units for the hidden layer for each of the obtained similarity groups. The blast data that were not used for training were used to validate the trained neural network models. For the same two similarity groups, multivariate regression models were also developed to predict mean particle size. Capability of the developed neural network models as well as multivariate regression models was determined by comparing predictions with measured mean particle size values and predictions based on one of the most applied fragmentation prediction models appearing in the blasting literature. Prediction capability of the trained neural network models as well as multivariate regression models was found to be strong and better than the existing most applied fragmentation prediction model. Diversity of the blasts data used is one of the most important aspects of the developed models. © Springer Science+Business Media B.V. 2012.
 Kulatilake, P. H., Wang, X., & Song, W. (2012). Stability investigations in threedimensions around a tunnel in a metal mine in China. 2012 SME Annual Meeting and Exhibit 2012, SME 2012, Meeting Preprints, 640652.More infoAbstract: To exploit an underground mine effectively and safely, it is important to have a good understanding of the geotechnical behavior around underground excavations made in the mine. The aim of this study was to investigate the geotechnical behavior around a tunnel excavated in a metal mine in China at the threedimensional level. Using the information available on lithology, geological structures, insitu stress measurements, physical and mechanical properties of intact rock, discontinuities and interfaces between different rocks, a threedimensional numerical model was built by using the 3DEC software to simulate a tunnel excavation under a high insitu stress condition. Effect of the discontinuity network, possible intact rock and discontinuity parameter variability, representation of rock masses as discontinuum or continuum material and rock support system on the deformation, stress, failure zone and stability around the tunnel was investigated. A very good agreement was obtained between the collected field deformation monitoring results and the results of the conducted numerical stress analyses. Copyright © 2012 by SME.
 Qiong, W. u., & Kulatilake, P. H. (2012). Application of equivalent continuum and discontinuum stress analyses in threedimensions to investigate stability of a rock tunnel in a dam site in China. Computers and Geotechnics, 46, 4868.More infoAbstract: Deformation and stability around a tunnel in a limestone rock mass is investigated using both (a) an equivalent continuum/discontinuum model and (b) a fully discontinuum model through threedimensional stress analyses. In model (a), the Representative Elementary Volume (REV) mechanical properties estimated in a previous paper are used to represent the combined behavior of intact rock and 4 of the five fracture sets that had significantly smaller sizes compared to the fifth fracture set (bedding planes). To this system the bedding plane fractures were added according to its statistical geometric properties. In model (b), all five fracture sets were generated according to the obtained fracture geometry statistical models. Using model (a), effect of several factors, such as REV property values, bedding plane modeling method, the average horizontal stress to the vertical stress and tunnel depth, on the geotechnical behavior around the tunnel was studied in detail. A comparison is made between the stress analyses results obtained through the aforementioned two models. © 2012 Elsevier Ltd.
 Qiong, W. u., & Kulatilake, P. H. (2012). REV and its properties on fracture system and mechanical properties, and an orthotropic constitutive model for a jointed rock mass in a dam site in China. Computers and Geotechnics, 43, 124142.More infoAbstract: Fracture data available for one of the rock masses (limestone) in the dam site of Yujian River Reservoir were used to build and validate a stochastic 3D fracture network model, and to perform a REV and equivalent continuum study in 3D. A number of relations are developed in the paper between the rock mass mechanical parameters and fracture tensor components in 3D. Based on the mechanical parameter values obtained in every 45° direction in 3D, the principal parameter values, principal directions and tensors are developed for rock mass mechanical parameters to represent the REV block size properties. An incrementally linear elastic, orthotropic constitutive model is suggested to represent the equivalent continuum prefailure mechanical behavior of the jointed rock mass by incorporating the effect of joint geometry network by the fracture tensor components. © 2012 Elsevier Ltd.
 Wang, L., Kulatilake, P. H., Tang, H., & Liang, Y. (2012). Rock slope stability study for Yujian River dam site based on kinematic analyses. Advanced Materials Research, 446449, 20482055.More infoAbstract: Lithological information, rock mass fracture data and discontinuity shear strength obtained through field investigations have been used in conducting kinematic analyses for the rock slopes that exist in the Yujian River dam site to evaluate the stability of the slopes. Results given in the paper can be considered as conservative because of several conservative assumptions used in the analyses. Dam site slopes seem quite stable up to 40° dip angle. Out of the three basic failure modes, possible wedge sliding seems to be the most likely one followed up with possible plane sliding as the second. Irrespective of the considered slope regions, slope dip direction ranges 270315° and 200210° seem to be the worst cases for possible instability of slopes in the dam site. Regional slopes in the dam site can be ranked with respect to safety from the lowest to highest in the following order: Rc1, Re2, Rc2, Rd1, Rb, Ra, Rd2 and Re1. Note that the dam site slopes are currently stable and the existing slope angles agree well with the results obtained from the rock slope stability analyses. © (2012) Trans Tech Publications, Switzerland.
 Wang, X., Kulatilake, P. H., & Song, W. (2012). Stability investigations around a mine tunnel through threedimensional discontinuum and continuum stress analyses. Tunnelling and Underground Space Technology, 32, 98112.More infoAbstract: To exploit an underground mine effectively and safely, it is important to have a good understanding of the geotechnical behavior around the underground excavations made in the mine. The aim of this study was to investigate the geotechnical behavior around a tunnel excavated in a metal mine in China at the threedimensional level. Using the information available on lithology, geological structures, in situ stress measurements, physical and mechanical properties of intact rock, discontinuities and interfaces between different rocks, a threedimensional numerical model was built by using the 3DEC software package to simulate a tunnel excavation under a high in situ stress condition. Effect of the discontinuity network, possible intact rock and discontinuity parameter variability, representation of rock masses as discontinuum or equivalent continuum material and rock support system on the deformation, and stability around the tunnel was investigated. A very good agreement was obtained between the collected field deformation monitoring results and the results of the conducted numerical stress analyses. © 2012 Elsevier Ltd.
 Zhang, Z. X., Xu, Y., Kulatilake, P. H., & Huang, X. (2012). Physical model test and numerical analysis on the behavior of stratified rock masses during underground excavation. International Journal of Rock Mechanics and Mining Sciences, 49, 134147.
 Zhengxing, Y. u., Kulatilake, P. H., & Jiang, F. (2012). Effect of Tunnel Shape and Support System on Stability of a Tunnel in a Deep Coal Mine in China. Geotechnical and Geological Engineering, 30(2), 383394.More infoAbstract: Stability level of tunnels that exist in an underground mine has a great influence on the safety, production and economic performance of the mine. Ensuring of stability for softrock tunnels is an important task for deep coal mines located in high in situ stress conditions. The aim of this study is to investigate the effect of tunnel shape and support pattern on the deformation, failure zone and stability around a tunnel located in a coal rock mass in China and to select an appropriate tunnel shape and a support pattern to provide a stable stressdeformation condition around the tunnel. Using the available information on stratigraphy, geological structures, in situ stress measurements and geomechanical properties of intact rock and discontinuity interfaces, a threedimensional numerical model was built using the FLAC software to simulate the stress conditions around the tunnel in the coal rock mass. Analyses were conducted for several tunnel shapes and rock support patterns. Results obtained for the distribution of failed zones, and stress and displacement fields around the tunnel were compared to select the best tunnel shape and support pattern to achieve the optimum stability conditions. Also, a comparison is given between the numerical predictions and field deformation monitoring results. © 2011 Springer Science+Business Media B.V.
 Hudaverdi, T., Kulatilake, P. H., & Kuzu, C. (2011). A new model for blast fragmentation prediction based on multivariate analysis. SME Annual Meeting and Exhibit and CMA 113th National Western Mining Conference 2011, 1115.
 Hudaverdi, T., Kulatilake, P. H., & Kuzu, C. (2011). Prediction of blast fragmentation using multivariate analysis procedures. International Journal for Numerical and Analytical Methods in Geomechanics, 35(12), 13181333.More infoAbstract: An extensive multivariate analysis procedure for prediction of blast fragmentation distribution is presented. Several blasts performed in various mines and rock formations in the world are brought together and evaluated. Blast design parameters, the modulus of elasticity, in situ block size are considered to perform multivariate analysis. The hierarchical cluster analysis is used to separate the blasts data into different groups of similarity. Group memberships were checked by the discriminant analysis. The multivariate regression analysis was applied to develop prediction equations for the estimation of the mean particle size of muckpiles. Two different prediction equations were developed based on the rock stiffness. Validation of the proposed equations on various mines is presented and the capability of the prediction equations was compared with one of the most applied fragmentation distribution models appearing in the blasting literature. Prediction capability of the proposed models was found to be strong. Diversity of the blasts data used is one of the most important aspects of the developed models. The models are not complex and suitable for practical use at mines.b © 2010 John Wiley & Sons, Ltd..
 Kulatilake, P. H., Wang, L. Q., Tang, H. M., & Liang, Y. (2011). Evaluation of rock slope stability for Yujian River dam site China by block theory analyses. 45th US Rock Mechanics / Geomechanics Symposium.More infoAbstract: Lithological information, rock mass fracture data and discontinuity shear strength obtained through field investigations have been used to conduct block theory analyses for the rock slopes that exist in the dam site to evaluate the stability of the slopes. The analyses were performed using mean discontinuity set orientations for each rock mass region under gravitational loading to calculate the maximum safe slope angles (MSSA) for different cut slope directions. The following conclusions have been made based on the block theory analysis results: (1) The final MSSA range between 30° and 47°, 44° and 70°, 47° and 69°for cut slope dip directions of 20°30°, 105°210°, and 270°355°, respectively; (2) For cut slope dip directions of 105°210° and 270°355°, wide ranges of values have been obtained for the final MSSA reflecting the influence of variability of fracture orientations on MSSA; (3) Apart from the region Rd1 for slope dip directions in the range 20°30°, rest of the regions at the dam site seem to be stable for slope angles less than 40°. © 2011 ARMA, American Rock Mechanics Association.
 Kulatilake, P. H., Wang, L., Tang, H., & Liang, Y. (2011). Evaluation of rock slope stability for yujian river dam site by kinematic and block theory analyses. Computers and Geotechnics, 38(6), 846860.More infoAbstract: Lithological information, rock mass fracture data and discontinuity shear strength obtained through field investigations have been used to conduct kinematic and block theory analyses for the rock slopes that exist in the dam site to evaluate the stability of the slopes. The analyses were performed using mean discontinuity set orientations for each rock mass region under gravitational loading to calculate the maximum safe slope angles (MSSA) for different cut slope directions. Results show that final MSSAs obtained from kinematic analysis are less than or equal to that obtained from block theory analysis. The following conclusions have been made based on the block theory analysis results, which are closer to the reality: (1) The final MSSA range between 30° and 47°, 44° and 70°, 47° and 69° for cut slope dip directions of 2030°, 105210°, and 270355°, respectively; (2) For cut slope dip directions of 2030°, 200210° and 275315°, wide ranges of values have been obtained for the final MSSA reflecting the influence of variability of fracture orientations on MSSA; (3) Apart from the region Rd1 for slope dip directions in the range 2030°, rest of the regions at the dam site seem to be stable for slope angles less than 40°. Detailed comparisons are given between the kinematic and block theory analyses covering both the theoretical concepts and application results. Also a brief comparison is included between the laboratory and in situ discontinuity shear strength results. © 2011 Elsevier Ltd.
 Nengxiong, X. u., Tian, H., Kulatilake, P. H., & Duan, Q. (2011). Building a three dimensional sealed geological model to use in numerical stress analysis software: A case study for a dam site. Computers and Geotechnics, 38(8), 10221030.More infoAbstract: The dam area of the SUOXI hydropower project shows high terrain undulation and complex geological conditions, containing 6 faults and 7weak interbeds. A geometric model developed to represent the geology and engineering structures should incorporate the geological realities and should allow suitable mesh generation to perform numerical stress analysis. This is an important precondition to perform rock mass stability analysis of a dam foundation based on a numerical stress analysis software such as FLAC3D. Using the modeling tools available in FLAC3D, it is difficult to construct a complex geological model even after performing a large amount of plotting and data analyses. The 3D geological modeling technique suggested in this paper, named as Sealed Geological Modeling (SGM), is a powerful tool for constructing complex geological models for rock engineering projects that require numerical stress analysis. Applying this technique, first, the geological interfaces are constructed for the dam area of SUOXI hydropower project using various interpolation procedures including geostatistical techniques. Then a unitary wire frame is constructed and the interfaces are connected seamlessly. As the next step, a block tracing technique is used to build a geological model that consists of 130 seamlessly connected blocks. Finally, based on the Advancing Front Technique (AFT), each block is discretized into tetrahedrons and a mesh is generated including 57,661 nodes and 215,471 tetrahedrons which is suitable to perform numerical stress analysis using FLAC3D. © 2011 Elsevier Ltd.
 Qiong, W. u., Kulatilake, P. H., & Tang, H. (2011). Comparison of rock discontinuity mean trace length and density estimation methods using discontinuity data from an outcrop in Wenchuan area, China. Computers and Geotechnics, 38(2), 258268.More infoAbstract: The equations that exist in the literature to estimate corrected mean trace length and corrected twodimensional density of a rock discontinuity set using area sampling technique are critically reviewed. The discontinuity traces appearing in an outcrop in Yingxiu area in China are used along with rectangular windows to calculate the corrected mean trace length and twodimensional density using Kulatilake and Wu's equations. Similarly, circular windows are used along with Mauldon's and Zhang and Einstein's equation to calculate the mean trace length and Mauldon's equation to calculate the twodimensional density for the same discontinuity sets using the same outcrop discontinuity trace data. For both parameters, the predictions based on the rectangular window methods were found to be more accurate than that based on the circular window methods. © 2010 Elsevier Ltd.
 Wang, L., Kulatilake, P. H., Tang, H., Liang, Y., & Qiong, W. u. (2011). Kinematic analyses of sliding and toppling failure of double free face rock mass slopes. Yantu Lixue/Rock and Soil Mechanics, 32(SUPPL. 1), 7277.More infoAbstract: Free face geometry shape plays an important role in determining failure and evaluating slope stability. Kinematic analysis is a very effective method in solving these issues stated above. Most kinematic analyses for rock mass slope stability mainly focus on the research of failure mode and maximum safe slope angle for single free face slopes. However, there are many slopes with double free face along strike of cut slope in engineering applications. In this paper, research on rock slope stability is extended to slopes with double free face. The corresponding sliding modes are divided into four types: (a) plane sliding along a discontinuity; (b) wedge sliding along a single discontinuity plane; (c) wedge sliding along the intersection line of two discontinuities and (d) toppling failure. On the stereographic projection, sliding zone for plane and wedge sliding can be marked as the area formed by the true dip direction lines of the two faces that form the double free face and the friction circle where the dip vector or the line of intersection vector should lie. The toppling zone can be marked as the area formed by the true dip direction lines of the two faces that form the double free face, friction circle and the periphery of the stereographic plot where the normal vector to the discontinuity should lie. The procedure to obtain the maximum safe slope angle for a double free face rock mass slope and the design principle to estimate the cut slope angle are suggested. Based on the method presented above, Shengjipo double free face rock mass slope, which is located in Zigui county, Hubei Province, is used as a case study to illustrate determination of failure modes and maximum safe slope angles.
 Wu, Q., Kulatilake, P. H., & Tang, H. M. (2011). Estimation of rock discontinuity mean trace length and density from an outcrop in Wenchuan, China. 45th US Rock Mechanics / Geomechanics Symposium.More infoAbstract: The equations exist in the literature to estimate corrected mean trace length and corrected two dimensional density of a discontinuity set using area sampling technique are critically reviewed. A shortcoming that exists with the equation proposed to calculate the corrected mean trace length of a discontinuity set based on a circular sampling window is pointed out. Also the major shortcomings that exist with the equation proposed to estimate the twodimensional density of a discontinuity set using a circular window are pointed out. The discontinuity traces appearing in the outcrop in Yingxiu area are used along with rectangular windows to calculate the corrected mean trace length and twodimensional density using Kulatilake and Wu's equations by incorporating the corrected relative frequency of each trace appearing on the two dimensional exposure. For both the corrected mean trace length and the twodimensional density, the predictions based on the rectangular window methods are found to be quite accurate. © 2011 ARMA, American Rock Mechanics Association.
 Wu, Q., Kulatilake, P. H., Hudaverdi, T., & Kuzu, C. (2011). A neural network approach to predict mean particle size in rock blast fragmentation. SME Annual Meeting and Exhibit and CMA 113th National Western Mining Conference 2011, 245254.More infoAbstract: Neural network methodology is used to predict mean particle size resulting from rock blast fragmentation. A blast data base developed in a previous study is used in the current study. A part of this blast data was used to train a singlehidden layer back propagation neural network model for each of the similarity groups obtained in the same previous study. LevenbergMarquardt algorithm provided the most stable and efficient training out of the four algorithms evaluated. An extensive analysis was performed using a part of the blast data to estimate the optimum number of units for the hidden layer for each similarity group. The remaining blast data are used to validate the trained neural network models. Capability of the developed neural network models is determined by comparing predictions with measured values and predictions based on one of the most applied fragmentation prediction models appearing in the blasting literature. Prediction capability of the neural network models was found to be strong and better than the existing most applied model. Diversity of the blasts data used is one of the most important aspects of the developed models. The developed neural network models are suitable for practical use at mines.
 Yu, Z. H., Kulatilake, P. H., & Jiang, F. X. (2011). Effect of tunnel shape and support system on stability of a tunnel in a deep coal mine in China. 45th US Rock Mechanics / Geomechanics Symposium.More infoAbstract: Stability level of tunnels that exist in an underground mine has a great influence on the safety, production and economic performance of the mine. Ensuring of stability for softrock tunnels is an important task for deep coal mines located in high in situ stress conditions. The aim of this study is to investigate the effect of tunnel shape and support pattern on the deformation, failure zone and stability around a tunnel located in a coal rock mass in China and to select an appropriate tunnel shape and a support pattern to provide stable stressdeformation conditions around the tunnel. Using the available information on stratigraphy, geological structures, in situ stress measurements and geomechanical properties of intact rock and discontinuity interfaces, a threedimensional numerical model was built by using FLAC software to simulate the stress conditions around the tunnel in the coal rock mass. Analyses were conducted for several tunnel shapes and rock support patterns. Results obtained for the distribution of failure zones, and stress and displacement fields around the tunnel were compared to select the best tunnel shape and support pattern to achieve the optimum stability conditions. © 2011 ARMA, American Rock Mechanics Association.
 Kulatilake, P. H. (2010). Preface: Special issue on Sri Lankan Geotechnical society's first international conference on soil & rock engineering. Geotechnical and Geological Engineering, 28(3), 209210.
 Kulatilake, P. H. (2010). Tensorial approach to rock mass strength and deformability in three dimensions. Analysis of Discontinuous Deformation: New Developments and Applications, 5971.
 Kulatilake, P. H., Qiong, W., Hudaverdi, T., & Kuzu, C. (2010). Mean particle size prediction in rock blast fragmentation using neural networks. Engineering Geology, 114(34), 298311.More infoAbstract: Multivariate analysis procedures and a neural network methodology are used to predict mean particle size resulting from rock blast fragmentation. A blast data base developed in a previous study is used in the current study. The data base consists of blast design parameters, explosive parameters, modulus of elasticity and insitu block size. In the same previous study a hierarchical cluster analysis was used to separate the blast data into two different groups of similarity based on the intact rock stiffness. In the same study the group memberships were confirmed by the discriminant analysis. A part of this blast data was used in this study to train a singlehidden layer backpropagation neural network model to predict mean particle size resulting from blast fragmentation for each of the obtained similarity groups. The mean particle size was considered to be a function of seven independent parameters. Four learning algorithms were considered to train neural network models. LevenbergMarquardt algorithm turned out to be the best one providing the highest stability and maximum learning speed. An extensive analysis was performed to estimate the optimum value for the number of units for the hidden layer for each of the obtained similarity groups. The blast data that were not used for training are used to validate the trained neural network models. For the same two similarity groups, multivariate regression models were also developed to predict mean particle size. Capability of the developed neural network models as well as multivariate regression models is determined by comparing predictions with measured mean particle size values and predictions based on one of the most applied fragmentation prediction models appearing in the blasting literature. Prediction capability of the trained neural network models as well as multivariate regression models was found to be strong and better than the existing most applied fragmentation prediction model. Diversity of the blasts data used is one of the most important aspects of the developed models. The developed neural network models and multivariate regression analysis models are suitable for practical use at mines. © 2010 Elsevier B.V.
 Kulatilake, P. H., Park, J., Balasingam, P., & Morgan, R. (2008). Quantification of aperture and relations between aperture, normal stress and fluid flow for natural single rock fractures. Geotechnical and Geological Engineering, 26(3), 269281.More infoAbstract: Accurate quantification of rock fracture aperture is important in investigating hydromechanical properties of rock fractures. Liquefied wood's metal was used successfully to determine the spatial distribution of aperture with normal stress for natural single rock fractures. A modified 3D box counting method is developed and applied to quantify the spatial variation of rock fracture aperture with normal stress. New functional relations are developed for the following list: (a) Aperture fractal dimension versus effective normal stress; (b) Aperture fractal dimension versus mean aperture; (c) Fluid flow rate per unit hydraulic gradient per unit width versus mean aperture; (d) Fluid flow rate per unit hydraulic gradient per unit width versus aperture fractal dimension. The aperture fractal dimension was found to be a better parameter than mean aperture to correlate to fluid flow rate of natural single rock fractures. A highly refined variogram technique is used to investigate possible existence of aperture anisotropy. It was observed that the scale dependent fractal parameter, Kv, plays a more prominent role than the fractal dimension, Da1d, on determining the anisotropy pattern of aperture data. A combined factor that represents both Da1d and Kv, Da1d × Kv, is suggested to capture the aperture anisotropy. © Springer Science+Business Media B.V. 2007.
 Wang, M., & Kulatilake, P. H. (2008). Understanding of hydraulic properties from configurations of stochastically distributed fracture networks. Hydrological Processes, 22(8), 11251135.More infoAbstract: A systematic investigation of the effect of configurations of stochastically distributed fracture networks on hydraulic behaviour for fractured rock masses could provide either quantitative or qualitative correlation between the structural configuration of the fracture network and its corresponding hydraulic behaviour, and enhance our understanding of appropriate application of groundwater flow and contaminant transport modelling in fractured rock masses. In this study, the effect of block sizes, intersection angles of fracture sets, standard deviations of fracture orientation, and fracture densities on directional block hydraulic conductivity and representative elementary volume is systematically investigated in two dimensions by implementing a numerical discrete fracture fluid flow model and incorporating stochastically distributed fracture configurations. It is shown from this investigation that the configuration of a stochastically distributed fracture network has a significant quantitative or qualitative effect on the hydraulic behaviour of fractured rock masses. Compared with the deterministic fracture configurations that have been extensively dealt with in a previous study, this investigation is expected to be more practical and adequate, since fracture geometry parameters are inherently stochastically distributed in the field. Moreover, the methodology and approach presented in this study may be generally applied to any fracture system in investigating the hydraulic behaviours from configurations of the fracture system while establishing a 'bridge' from the discrete fracture network flow modelling to equivalent continuum modelling in fractured rock masses. Copyright © 2007 John Wiley & Sons, Ltd.
 Kulatilake, P. H., Park, J., Balasingam, P., & Mckenna, S. A. (2007). Hierarchical probabilistic regionalization of volcanism for Sengan region, Japan. Geotechnical and Geological Engineering, 25(1), 79102.More infoAbstract: A 1 km square regular grid system created on the Universal Transverse Mercator zone 54 projected coordinate system is used to work with volcanism related data for Sengan region. The following geologic variables were determined as the most important for identifying volcanism: geothermal gradient, groundwater temperature, heat discharge, groundwater pH value, presence of volcanic rocks and presence of hydrothermal alteration. Data available for each of these important geologic variables were used to perform directional variogram modeling and kriging to estimate geologic variable vectors at each of the 23949 centers of the chosen 1 km cell grid system. Cluster analysis was performed on the 23949 complete variable vectors to classify each center of 1 km cell into one of five different statistically homogeneous groups with respect to potential volcanism spanning from lowest possible volcanism to highest possible volcanism with increasing group number. A discriminant analysis incorporating Bayes' theorem was performed to construct maps showing the probability of group membership for each of the volcanism groups. The said maps showed good comparisons with the recorded locations of volcanism within the Sengan region. No volcanic data were found to exist in the group 1 region. The high probability areas within group 1 have the chance of being the no volcanism region. Entropy of classification is calculated to assess the uncertainty of the allocation process of each 1 km cell center location based on the calculated probabilities. The recorded volcanism data are also plotted on the entropy map to examine the uncertainty level of the estimations at the locations where volcanism exists. The volcanic data cell locations that are in the high volcanism regions (groups 4 and 5) showed relatively low mapping estimation uncertainty. On the other hand, the volcanic data cell locations that are in the low volcanism region (group 2) showed relatively high mapping estimation uncertainty. The volcanic data cell locations that are in the medium volcanism region (group 3) showed relatively moderate mapping estimation uncertainty. Areas of high uncertainty provide locations where additional site characterization resources can be spent most effectively. The new data collected can be added to the existing database to perform future regionalized mapping and reduce the uncertainty level of the existing estimations. © Springer Science+Business Media B.V. 2006.
 Kulatilake, P. H., Balasingam, P., Park, J., & Morgan, R. (2006). Natural rock joint roughness quantification through fractal techniques. Geotechnical and Geological Engineering, 24(5), 11811202.More infoAbstract: Accurate quantification of roughness is important in modeling hydromechanical behavior of rock joints. A highly refined variogram technique was used to investigate possible existence of anisotropy in natural rock joint roughness. Investigated natural rock joints showed randomly varying roughness anisotropy with the direction. A scale dependant fractal parameter, Kv, seems to play a prominent role than the fractal dimension, Dr1d, with respect to quantification of roughness of natural rock joints. Because the roughness varies randomly, it is impossible to predict the roughness variation of rock joint surfaces from measurements made in only two perpendicular directions on a particular sample. The parameter D r1d × Kv seems to capture the overall roughness characteristics of natural rock joints well. The onedimensional modified divider technique was extended to two dimensions to quantify the twodimensional roughness of rock joints. The developed technique was validated by applying to a generated fractional Brownian surface with fractal dimension equal to 2.5. It was found that the calculated fractal parameters quantify the rock joint roughness well. A new technique is introduced to study the effect of scale on twodimensional roughness variability and anisotropy. The roughness anisotropy and variability reduced with increasing scale. © Springer 2006.
 Kulatilake, P. H., Park, J., & Malama, B. (2006). A new rock mass failure criterion for biaxial loading conditions. Geotechnical and Geological Engineering, 24(4), 871888.More infoAbstract: To simulate brittle rocks, a mixture of glastone, sand and water was used as a model material. Thin galvanized sheets of thickness 0.254 mm were used to create joints in blocks made out of the model material. To investigate the failure modes and strength, both the intact material blocks as well as jointed model material blocks of size 35.6 × 17.8 × 2.5 cm having different joint geometry configurations were subjected to uniaxial and biaxial compressive loadings. A new intact rock failure criterion is proposed at the 3D level. This criterion is validated for biaxial loading through laboratory experimental results obtained on intact model material blocks. Results obtained from both the intact and jointed model material blocks are used to develop a strongly nonlinear new rock mass failure criterion for biaxial loading. In this failure criterion, the fracture tensor component is used to incorporate the directional effect of fracture geometry system on jointed block strength. The failure criterion shows the important role, the intermediate principal stress plays on rock mass strength. © Springer 2006.
 Kulatilake, P. H., Park, J., Balasingam, P., & Morgan, R. (2006). Natural rock fracture aperture properties through fractals. Proceedings of the 41st U.S. Rock Mechanics Symposium  ARMA's Golden Rocks 2006  50 Years of Rock Mechanics.More infoAbstract: Accurate quantification of rock fracture aperture is important in investigating hydromechanical properties of rock fractures. Liquefied wood's metal was used successfully to determine the spatial distribution of aperture with normal stress for natural single rock fractures. A modified 3D box counting method is developed and applied to quantify the spatial variation of rock fracture aperture with normal stress. New functional relations are developed for the following list: (a) Aperture fractal dimension versus effective normal stress; (b) Aperture fractal dimension versus mean aperture; (c) Fluid flow rate per unit hydraulic gradient per unit width versus mean aperture; (d) Fluid flow rate per unit hydraulic gradient per unit width versus aperture fractal dimension. The aperture fractal dimension was found to be a better parameter than mean aperture to correlate to fluid flow rate of natural single rock fractures. A highly refined variogram technique is used to investigate possible existence of aperture anisotropy. It was observed that the scale dependent fractal parameter, Kv, plays a more prominent role than the fractal dimension, Dald, on determining the anisotropy pattern of aperture data. A combined factor that represents both Dald and Kv, Dald × Kv, is suggested to capture the aperture anisotropy. Copyright 2006, ARMA, American Rock Mechanics Association.
 Kulatilake, P. H., Park, J., & Um, J. (2004). Estimation of rock mass strength and deformability in 3D for a 30 m cube at a depth of 485 m at Äspö Hard Rock Laboratory. Geotechnical and Geological Engineering, 22(3), 313330.More infoAbstract: The rock fracture data provided by Swedish Nuclear Fuel and Waste Management Company were used to develop a 3D stochastic fracture network model for a 30m cube of Äspö diorite located at a depth of 485 m at Äspö Hard Rock Laboratory, Sweden. This fracture network model was validated. A new procedure is developed to estimate rock block strength and deformability in threedimensions allowing for the anisotropy and incorporating the inherently statistical fracture geometry for the selected cube. The mean rock mass strength was found to be 47% of the mean intact rock strength of 297 MPa at 485 m depth. The mean rock mass modulus was found to be 51% of the intact rock Young's modulus of 73 GPa. The rock mass Poisson's ratio was found to be 21% higher than the intact rock Poisson's ratio of 0.28. These percentages indicate the level of weakening of the rock mass due to the presence of fractures. The ratio of mean major principal rock mass strength/mean minor principal rock mass strength turned out to be 1.28. The ratio of mean major principal rock mass modulus/mean minor principal rock mass modulus turned out to be 1.21. © 2004 Kluwer Academic Publishers.
 Kulatilake, P. H., & Um, J. (2003). Spatial variation of cone tip resistance for the clay site at Texas A & M University. Geotechnical Special Publication, 4160.More infoAbstract: The cone tip resistance data available for the clay site at Texas A & M University (one of the National Geotechnical Experimentation Sites) are used to show how the spatial variation of a soil property can be quantified. It is suggested that the first step in quantifying the spatial variation of a soil property should be the identification or selection of the statistically homogeneous soil layers. A new simple procedure is suggested to identify statistically homogeneous layers in a soil profile. Through examples it is shown that the procedure works well in identifying the statistically homogeneous layers. For the chosen statistically homogeneous layers, the spatial variation of cone tip resistance in the depthwise is quantified in terms of a constant mean or a mean trend, variance/standard deviation/coefficient of variation or a standard deviation around the mean trend, and a correlation or variogram function. Correlation distances in the depth direction were ranged between 0.1 and 0.5m for the two soil layers investigated. It was shown that the correlation distance can decrease in the presence of a global mean trend for the soil property. In such cases, it is important to note that a part of the correlation is automatically included in the mean trend function. Cone penetrometer data available at spacings greater than or equal to 1m failed to indicate presence of correlation structures in the horizontal direction. Selfaffine fractals seem to have potential in capturing the different autocorrelation structures that exist with spatial variation of soil properties.
 Kulatilake, P. H., & Um, J. (2003). Spatial variation of cone tip resistance for the clay site at Texas A and M University. Geotechnical and Geological Engineering, 21(2), 149165.More infoAbstract: The cone tip resistance data available for the clay site at Texas A and M University, USA (one of the National Geotechnical Experimentation Sites) are used to show how the spatial variation of a soil property can be quantified. It is suggested that the first step in quantifying the spatial variation of a soil property should be the identification or selection of the statistically homogeneous soil layers. A new simple procedure is suggested to identify statistically homogeneous layers in a soil profile. Through examples it is shown that the procedure works extremely well in identifying the statistically homogeneous layers. For the chosen statistically homogeneous layers, the spatial variation of cone tip resistance with depth is quantified in terms of a constant mean or a mean trend, variance/standard deviation/coefficient of variation or a variance around the mean trend, and a correlation or variogram function. Correlation distances in the depth direction were found to be between 0.1 and 0.5 m for the two soil layers investigated. It was shown that the correlation distance decreases in the presence of a global mean trend for the soil property. In such cases, it is important to note that a part of the correlation is automatically included in the mean trend function.
 Kulatilake, P. H., Um, J., Wang, M., Escandon, R. F., & Narvaiz, J. (2003). Stochastic fracture geometry modeling in 3D including validations for a part of Arrowhead East Tunnel, California, USA. Engineering Geology, 70(12), 131155.More infoAbstract: Eighthundred and fifty nine fractures of a gneissic rock mass were mapped using 16 scanlines placed on steep rock exposures that were within 300 m of a tunnel alignment before the tunnel excavation. These data were analyzed using the software package FRACNTWK to find the number of fracture sets that exist in the rock mass, 3D fracture frequency for each set and the probability distributions of orientation, trace length, fracture size in three dimensions (3D) and spacing for each of the fracture sets. In obtaining these distributions corrections were applied for sampling biases associated with orientation, trace length, size and spacing. Developed stochastic 3D fracture network for the rock mass was validated by comparing statistical properties of observed fracture traces on the scanlines with the predicted fracture traces on similar scanlines. The onedimensional (1D) fracture frequency of the rock mass in all directions in 3D was calculated and is presented in terms of a stereographic plot. The 1D fracture frequency along the tunnel alignment direction was predicted to be about 6.5 fractures/m before the tunnel excavation. This prediction was found to be in excellent agreement with the observed values obtained about 1 year later during the tunnel excavation. This was another validation conducted for the developed 3D fracture network. © 2003 Elsevier Science B.V. All rights reserved.
 Malama, B., & Kulatilake, P. H. (2003). Models for normal fracture deformation under compressive loading. International Journal of Rock Mechanics and Mining Sciences, 40(6), 893901.More infoAbstract: A new semiempirical model that can be used to predict fracture deformation behavior under normal compressive loading is presented. The development of a simple exponential model is presented first after which a modified and more general exponential model, with an additional degree of freedom in the model parameters, is obtained. The simple and the modified exponential models are then compared to available fracture closure models, namely the empirical BartonBandis hyperbolic model, and a powerlaw model based on Hertzian contact theory, to determine how good they fit the results of fracture closure experiments under monotonically increasing normal compressive loading. A new parameter called the halfclosure stress, σ1/2, is introduced and is used, in addition to the maximum fracture closure, Δvm, in the model fitting procedures for the BartonBandis and the simple exponential model. The halfclosure stress is shown to be related to the initial normal stiffness, Kni, used in the original BartonBandis model. An additional parameter, n, is used in fitting the modified exponential model to the experimental data. Of the models presented herein, the modified exponential model was found to provide the best fit to the experimental data, for the same values of σ1/2 and Δvm, over the entire range of compressive stresses. The powerlaw model based on Hertzian contact theory was found to be unsuitable for predicting normal fracture deformation behavior. © 2003 Elsevier Ltd. All rights reserved.
 Zhang, Z., & Kulatilake, P. H. (2003). A new stereoanalytical method for determination of removal blocks in discontinuous rock masses. International Journal for Numerical and Analytical Methods in Geomechanics, 27(10), 791811.More infoAbstract: The paper provides a new stereoanalytical method, which is a combination of the stereographic method and analytical methods, to separate finite removable blocks from the infinite and tapered blocks in discontinuous rock masses. The methodology has applicability to both convex and concave blocks. Application of the methodology is illustrated through examples. Addition of this method to the existing block theory procedures available in the literature improves the capability of block theory in solving practical problems in rock engineering. © 2003 John Wiley and Sons, Ltd.
 Wang, M., Kulatilake, P. H., Um, J., & Narvaiz, J. (2002). Estimation of REV size and threedimensional hydraulic conductivity tensor for a fractured rock mass through a single well packer test and discrete fracture fluid flow modeling. International Journal of Rock Mechanics and Mining Sciences, 39(7), 887904.More infoAbstract: A new methodology is presented to determine the representative elementary volume (REV) size and threedimensional (3D) hydraulic conductivity tensor for a fractured rock mass. First, a 3D stochastic fracture network model was built and validated for a gneissic rock mass based on the fracture data mapped from scanline surveys at the site. This validated fracture network model was combined with the fracture data observed on a borehole to generate a stochasticdeterministic fracture network system in a cubic block around each packer test conducted at a different depth region in the same borehole. Each packer test was simulated numerically applying a developed discrete fracture fluid flow model to estimate the influenced region or effective range for the packer test. A cubic block of size 18 m, with the packer test interval of length about 6.5 m located at the centre of this block, was found to be suitable to represent the influenced region. Using this block size, the average flow rate per unit hydraulic gradient (defined as the transmissivity multiplied by mean width of flow paths) field for fractures was calibrated at different depth regions around the borehole by numerically simulating the packer tests conducted at different depth regions. The average flow rate per unit hydraulic gradient of the fractures that intersect the borehole was considered to be quite different to the average flow rate per unit hydraulic gradient of the fractures that do not intersect the borehole. A relation was developed to quantity the ratio between these two parameters. By studying the directional hydraulic conductivity behaviour of different cubic block sizes having the validated stochastic fracture network and calibrated hydraulic parameters, a REV for the hydraulic behaviour of the rock mass was estimated to be a block size of 15 m. The hydraulic conductivity tensor in 3D computed through regression analysis using the calculated directional hydraulic conductivity values in many directions was found to be significantly anisotropic. The principal directions of the hydraulic conductivity tensor were found to be agreeable with the existing fracture system of the site. Further, the geometric hydraulic conductivity calculated was found to be comparable to the hydraulic conductivity estimated through the radial flow assumption in continuum porous media. © 2002 Elsevier Science Ltd. All rights reserved.
 Kulatilake, P. H., Liang, J., & Gao, H. (2001). Experimental and numerical simulations of jointed rock block strength under uniaxial loading. Journal of Engineering Mechanics, 127(12), 12401247.More infoAbstract: To simulate brittle rocks, a mixture of sand, plaster of paris, and water was used as a model material. Thin galvanized sheets were used to create joints in blocks made out of the model material. To investigate the failure modes and strength, 30 × 12.5 × 8.6 cm jointed model material blocks having different joint geometry configurations were subjected to uniaxial compressive loading. Results indicated three failure modes: (1) tensile failure through intact material; (2) combined shear and tensile failure or only shear failure on joints; and (3) mixed failure of the above two modes depending on the joint geometry. The fracture tensor component in a certain direction quantifies the directional effect of the joint geometry, including number of fracture sets, fracture density, and probability distributions for size and orientation of these fracture sets. Results obtained from the experiments were used to develop a strongly nonlinear relation between the fracture tensor component and the jointed block strength. The laboratory experiments conducted on jointed model material blocks were simulated numerically using the Universal Distinct Element Code (UDEC). With careful selection of suitable material constitutive models for intact model material and model joints, and accurate estimation and calibration of mechanical parameters of the constitutive models through a combination of laboratory testing and numerical simulations of the intact model material and model joints separately, it was possible to obtain a good agreement between the laboratory experimental and distinct element numerical results.
 Kulatilake, P. H., Malama, B., & Wang, J. (2001). Physical and particle flow modeling of jointed rock block behavior under uniaxial loading. International Journal of Rock Mechanics and Mining Sciences, 38(5), 641657.More infoAbstract: Laboratory experiments and numerical simulations, using Particle Flow Code (PFC3D), were performed to study the behavior of jointed blocks of model material under uniaxial loading. The effect of joint geometry parameters on the uniaxial compressive strength of jointed blocks was investigated and this paper presents the results of the experiments and numerical simulations. The fracture tensor component in a given direction is used to quantify the combined directional effect of joint geometry parameters including joint density, orientation and size distributions, and the number of joint sets. The variation of the uniaxial compressive strength of the jointed blocks of the model material with the fracture tensor component was investigated. Both the laboratory experiments and the numerical simulations showed that the uniaxial block strength decreases in a nonlinear manner with increasing values of the fracture tensor component. It was observed that joint geometry configuration controls the mode of failure of the jointed blocks and three modes of failure were identified, namely (a) tensile splitting through the intact material, (b) failure by sliding along the joint plane and/or by displacement normal to the joint plane and, (c) mixed mode failure involving both the failure mechanisms in (a) and (b). It has also been shown that with careful parameter calibration procedures, PFC3D could be used to model the strength behavior of jointed blocks of rock. © 2001 Elsevier Science Ltd. All rights reserved.
 Um, J. ., & Kulatilake, P. H. (2001). Kinematic and block theory analyses for shiplock slopes of the three Gorges dam site in China. Geotechnical and Geological Engineering, 19(1), 2142.More infoAbstract: The granitic rock mass that exists in the shiplock region of the Three Gorges dam site contains a number of major discontinuities and about four sets of minor discontinuities. One hundred and thirty three major discontinuities have been mapped around the shiplock covering an area of 1740 × 600 m. These major discontinuities were used to perform rock slope kinematic and block theory analyses. Kinematic analyses were performed under the following two cases: (1) assuming all the mapped discontinuities cross the shiplock; (2) using only discontinuities that actually intersect the shiplock. Under case (1) and case (2) the shiplock faces in the proposed permanent shiplock region in fresh rock were found to be stable up to a cut slope of about 45° and 58°, respectively. Block theory was applied to identify different block types that exist on the shiplock faces and to estimate the maximum safe slope angles on the shiplock faces. The orientations of the major discontinuities that actually intersect the shiplocks were considered in this analysis. The total length of the shiplock (1750 m) was divided into 50 m segments. From the stereoplots, the key blocks (Type I) and/or potential key blocks (Type I) were found for only five segments of the shiplock slopes. It was found that the dip of the cut slope should be less than about 60° to avoid creation of a key block on the proposed shiplock slopes. However, it is important to keep in mind that these conclusions are based on the kinematic analyses performed using only the major discontinuities. Further kinematic as well as kinetic analyses are recommended incorporating minor discontinuities, water forces, earthquake forces etc. before making the final conclusions about the maximum safe slope angle for the shiplock region.
 Wang, M., Kulatilake, P. H., Panda, B. B., & Rucker, M. L. (2001). Groundwater resources evaluation case study via discrete fracture flow modeling. Engineering Geology, 62(4), 267291.More infoAbstract: The Arizona Department of Transportation (ADOT) is preparing to upgrade State Route 260 between Payson and Heber. It is estimated that a total of about one million cubic meters of water will be required for embankment construction during a period of about 84 months to upgrade the first 33.8 km of the highway. ADOT is investigating various sources of construction water for use in the highway improvement project, including groundwater resources along the highway corridor. A region known as the RV site, underlain by fractured granite, is located 12.9 km east of Payson. The site includes three springs, a creek and several wells. Several boreholes and observation wells were made to a maximum depth of 157 m to obtain fracture data and to conduct pumping tests with monitoring. Fracture data recorded by acoustic televiewer logs were used to build a fracture network model for the rock mass. Results of a 24hour and a 7day pumping tests were used to calibrate hydraulic parameters of a finite element discrete fracture fluid flow model considering the region as a heterogeneous, anisotropic, fractured medium. A 38day multiwell pumping test was used to validate the calibrated numerical model. The calibrated model showed the capability to provide reasonably accurate predictions for new pumping tests conducted in the same well field. The validated model was then used to simulate pumping exceeding a 7 year period under different scenarios incorporating different sets of boundary conditions and different pumping rates at multiwells, with and without recharge, to evaluate the yield of the aquifer and to assess the effect of longterm pumping on the environment. The results indicated that (a) the combined yield of the wells in the RV site is sufficient to meet the water demand for the ADOT highway project and (b) the water levels in the well field would decline between 3.0 and 7.6 m after one year of pumping and by 12.2 to over 30.5 m during the life of the project. © 2001 Elsevier Science B.V. All rights reserved.
 Kulatilake, P. H., & Panda, B. B. (2000). Effect of block size and joint geometry on jointed rock hydraulics and REV. Journal of Engineering Mechanics, 126(8), 850858.More infoAbstract: The effect of joint geometry parameters and the size of jointed rock block on the hydraulic properties of jointed rock, including the equivalent continuum behavior, was investigated through numerical experimentation. The chance to reach equivalent continuum behavior for a rock mass having a certain joint configuration increases with the increase of block size. The REV (representative elementary volume) size to show hydraulic equivalent continuum behavior for a jointed rock system was found to depend on the orientation of the joint sets; the REV size seems to decrease with increasing joint density and joint size. The REV does not seem to exist for some rock masses having joint systems with low relative orientation angles (systems with two joint sets) and low densities. The average block permeability (K0) value at the REV size for a joint system increases with increase in joint size and joint density. The equivalent continuum behavior of a joint system can be expressed with respect to a cuttoff value for the first invariant of fracture tensor (F0). The block size corresponding to the aforementioned cutoff F0 (between 10 and 30 for the joint systems investigated in this study) can be considered to provide the REV size for a given joint configuration. A threshold value for F0 (a joint network having two perpendicular joint sets produced 2.75) can be used to find the chance for a given joint network to have nonzero block permeability. A strong power functional relation seems to exist between the directional permeability and the fracture tensor component for the connected joint configuration when rock blocks contain minor discontinuities.
 Kulatilake, P. H., & Um, J. (1999). Requirements for accurate quantification of selfaffine roughness using the roughnesslength method. International Journal of Rock Mechanics and Mining Sciences, 36(1), 518.More infoAbstract: Selfaffine fractals have the potential to represent rock joint roughness profiles. Fractional Brownian profiles (selfaffine profiles) with known values of fractal dimension, D, input standard deviation, σ, and data density, d, were generated. For different values of the input parameter of the roughnesslength method (window length, w), D and another associated fractal parameter A were calculated for the aforementioned profiles. The calculated D was compared with the D used for the generation to determine the accuracy of calculated D. Suitable ranges for w were estimated to produce accurate D (within ±10% error) for the generated profiles. The results showed that to obtain reliable estimates for fractal parameters of a natural rock joint profile, it is necessary to choose a unit for the profile length to satisfy a data density (d) greater than or equal to 5.1. For roughness profiles having 5.1 ≤ d ≤ 51.23 and 1.2 ≤ D ≤ 1.7, w values between 2.5% and 10% of the profile length were found to be highly suitable to produce accurate fractal parameter estimates. It is recommended to use at least seven w values between the estimated minimum and maximum suitable w values in estimating fractal parameters of a natural rock joint profile. It was found that σ and a global trend of a roughness profile have no effect on calculated D. The estimated A was found to increase with both D and σ. The parameter D captures the autocorrelation and A captures the amplitude of a roughness profile at different scales. Therefore, the parameters D and A are recommended to use with the roughnesslength method in quantifying rock joint roughness. In addition, at least one more parameter is required to quantify the global trend of a roughness profile, if it exist; in many cases just the inclination or declination angle of the roughness profile in the direction considered would be sufficient to estimate the global trend. Calculated crossover lengths (segment length of a profile at which a selfaffine profile becomes selfsimilar) for the profiles investigated were found to be extremely small (less than 0.6% of the profile length) indicating that laser profilometers are required to make roughness measurements at interval lengths comparable to the crossover lengths of the natural rock joint profiles. To calculate rock joint roughness parameters accurately using the selfsimilar techniques, it is necessary to have roughness measurements made at interval lengths comparable to the crossover length of the profile. This indicate clearly the difficulty of using selfsimilar techniques such as the divider method in estimating rock joint roughness accurately.
 Kulatilake, P. H., Um, J., Panda, B. B., & Nghiem, N. (1999). Development of new peak shearstrength criterion for anisotropic rock joints. Journal of Engineering Mechanics, 125(9), 10101017.More infoAbstract: The roughness of a natural rock joint was measured in different directions using a laser profilometer. Two stationary roughness parameters and a nonstationary roughness parameter (all fractal based) were used to quantify anisotropic roughness. A plaster of Paris based model material was used to make model material replicas of the natural rock joint. Direct shear tests were performed at five different normal stresses, in each of the directions that were used for the roughness measurements, to measure the anisotropic peak shear strength of the model joint. Required observations and experiments were conducted to estimate (1) the asperity shear area as a proportion of the total surface area of the joint, for each tested joint; (2) the basic friction angle of the model material; and (3) the joint compressive strength. Tests were also conducted to develop a peak shearstrength criterion for the intact model material. Part of the direct shear test data was used to develop a new peak shearstrength criterion for joints including the aforementioned parameters. The other part of the data was used for model validation.
 Panda, B. B., & Kulatilake, P. H. (1999). Effect of joint geometry and transmissivity on jointed rock hydraulics. Journal of Engineering Mechanics, 125(1), 4150.More infoAbstract: The effect of joint geometry parameters (joint orientation, joint density, and joint size) and joint transmissivity on the hydraulic properties of jointed rock blocks was examined in two dimensions using joint networks with two joint sets through numerical experiments. The chance for equivalent continuum hydraulic behavior and the average block permeability of jointed rock were found to be very sensitive to the relative orientation of the joint sets in a joint network system. Although both the equivalent continuum behavior and the average block permeability seem to be little affected by the change of the standard deviation of the joint orientation distribution, the degree of permeability anisotropy seems to depend on the standard deviation of the joint orientation especially for anisotropic joint systems. Chance for the equivalent continuum behavior and the average block permeability were found to increase with the increase of joint density and joint size. The chance for equivalent continuum behavior and the permeability anisotropy were found to depend on the distribution of joint hydraulic conductivity. The influence of the joint hydraulic conductivity is more prominent for joint systems with high joint density and large joint sizes. The greater the variation of the hydraulic conductivity value among the joints in a joint set, the lesser the chance for the medium to behave like an equivalent continuum porous medium.
 Panda, B. B., & Kulatilake, P. H. (1999). Relations between fracture tensor parameters and jointed rock hydraulics. Journal of Engineering Mechanics, 125(1), 5159.More infoAbstract: The fracture tensor is a combined measure of joint geometry parameters. The possible relations between the properties of the fracture tensor and the hydraulic properties of jointed rock including the equivalent continuum behavior are investigated in this paper through numerical experiments performed at the twodimensional level incorporating joint networks with two joint sets. The possibility of a nonzero average block permeability and the equivalent continuum behavior of a joint system can be expressed with the help of the first invariant of fracture tensor (F0). A threshold value of F0 (to achieve nonzero block permeability) between 3 and 6 and a cutoff value of F0 (to achieve equivalent continuum behavior) between 10 and 30 were obtained for the joint geometry configurations studied in the investigation. Both the threshold value and the cutoff value of F0 were found to depend on the relative orientation of joint sets. Joint systems with low relative orientation angles (for systems with only two joint sets) and low densities do not show hydraulically equivalent continuum behavior. A nonlinear relation seems to exist between the directional permeability and the fracture tensor component for the connected joint configuration when rock blocks contain only minor discontinuities. For rock blocks containing major discontinuities, this relation becomes linear.
 Kulatilake, P. H., Um, J., & Pan, G. (1998). Requirements for accurate quantification of selfaffine roughness using the variogram method. International Journal of Solids and Structures, 35(3132), 41674189.More infoAbstract: Both stationary and nonstationary fractional Brownian profiles (selfaffine profiles) with known values of fractal dimension, D, input standard deviation, σ, and data density, d, were generated. For different values of the input parameter of the variogram method (lag distance, h), D and another associated fractal parameter Kv were calculated for the aforementioned profiles. It was found that σ has no effect on calculated D. The estimated Kv was found to increase with D, σ and d according to the equation Kv = 2.0 × 105d0.35σ1.95D14.5. The parameter Kv seems to have potential to capture the scale effect of roughness profiles. Suitable ranges for h were estimated to obtain computed D within ±10% of the D used for the generation and also to satisfy a power functional relation between the variograrn and h. Results indicated the importance of removal of nonstationarity of profiles to obtain accurate estimates for the fractal parameters. It was found that at least two parameters are required to quantify stationary roughness; the parameters D and Kv are suggested for use with the variograrn method. In addition, one or more parameters should be used to quantify the nonstationary part of roughness, if it exists; at the basic level, the mean inclination/declination angle of the surface in the direction considered can be used to represent the nonstationarity. A new concept of feature size range of a roughness profile is introduced in this paper. The feature size range depends on d, D and σ. The suitable h range to use with the variogram method to produce accurate fractal parameter values for a roughness profile was found to depend on both d and D. It is shown that the feature size range of a roughness profile plays an important role in obtaining accurate roughness parameter values with both the divider and the variogram methods. The minimum suitable h was found to increase with decreasing d and increasing D. Procedures are given in this paper to estimate a suitable h range for a given natural rock joint profile to use with the variogram method to estimate D and Kv accurately for the profile. © 1998 Elsevier Science Ltd. All rights reserved.
 H., P., Fiedler, R., & Panda, B. B. (1997). Box fractal dimension as a measure of statistical homogeneity of jointed rock masses. Engineering Geology, 48(34), 217229.More infoAbstract: A written computer programme to estimate the box fractal dimension (DB) is verified by estimating DB of the triadic Koch curve for which the theoretical D is known. The influence of a number of input parameters of the boxcounting method on the accuracy of estimated DB is evaluated using the same Koch curve. The employed size range of the applied box networks was found to be the parameter which has the strongest influence on the accuracy of estimated DB. This indicated the importance of finding the range of selfsimilarity or selfaffinity for the object considered to select the proper range for the box sizes and, in turn, to obtain accurate estimates of DB. By calculating DB for different block sizes sampled from three generated twodimensional joint patterns, it is shown that DB can capture the combined effect of jointsize distribution and joint density on the statistical homogeneity of rock masses. The spatial variation of DB along a 350 m stretch of a tunnel in the shiplock area of the Three Gorges dam site is computed using the joint data mapped on the walls and the roof of the tunnel. This spatial variation of DB is used, along with the visual geological evaluation of the joint trace maps of the tunnel, in making decisions about statistical homogeneity of the rock mass around the tunnel. The results obtained on statistically homogeneous regions were found to be quite similar to the results obtained from a previous statistical homogeneity investigation which incorporated the effect of number of joint sets and their orientation distributions, but not the spatial variation of DB. It is recommended that the spatial variation of DB is used, along with the results of other methods such as contingency table analysis and equal area plots, which incorporate the effect of joint orientation distribution, in addition to the geology of the site, in determining the statistically homogeneous regions of jointed rock masses. © 1997 Elsevier Science B.V.
 Kulatilake, P. H., & Um, J. (1997). Requirements for accurate quantification of self affine roughness using the roughnesslength method. International journal of rock mechanics and mining sciences & geomechanics abstracts, 34(34), 392393.More infoAbstract: Accurate quantification of roughness is important in modeling strength, deformability and fluid flow behaviors of joints. Selfaffine fractals seem to have potential to represent rock joint roughness profiles. Both stationary and nonstationary fractional Brownian profiles with known values of fractal dimension, D, and input standard deviation, σ, were generated. For different values of the input parameter of the roughnesslength method, D and another associated fractal parameter A were calculated for the aforementioned profiles. In general, significant results were obtained.
 Kulatilake, P. H., He, W., Um, J., & Wang, H. (1997). A physical model study of jointed rock mass strength under uniaxial compressive loading. International journal of rock mechanics and mining sciences & geomechanics abstracts, 34(34), 692693.More infoAbstract: Jointed rock mass strength is investigated through physical modeling using model material blocks having different joint configurations subjected to only uniaxial compressive loading. Jointed model material blocks exhibited different failure modes depending on the joint configurations. Orientation of the joint sets played a significant role related to the modes of failure (Fig. 1). It was possible to obtain a strong nonlinear relation between the jointed model mass strength and the fracture tensor component to cover the strengths resulting from all the different failure modes observed in the investigation (Fig. 2). The fracture tensor component was used to obtain the combined effect of a number of joint sets, joint density, and distributions of size and orientation of the joint sets in a chosen direction. Future research is suggested to generalize the promising rock mass strength criterion obtained in the performed research.
 Kulatilake, P. H., Um, J., & Pan, G. (1997). Requirements for accurate estimation of fractal parameters for selfaffine roughness profiles using the line scaling method. Rock Mechanics and Rock Engineering, 30(4), 181206.More infoAbstract: A new concept of feature size range of a roughness profile is introduced in the paper. It is shown that this feature size range plays an important role in estimating the fractal dimension, D, accurately using the divider method. Discussions are given to indicate the difficulty of using both the divider and the box methods in estimating D accurately for selfaffine profiles. The line scaling method's capability in quantifying roughness of natural rock joint profiles, which may be selfaffine, is explored. Fractional Brownian profiles (selfaffine profiles) with and without global trends were generated using known values of D, input standard deviation, σ, and global trend angles. For different values of the input parameter of the line scaling method (step size a0), D and another associated fractal parameter C were calculated for the aforementioned profiles. Suitable ranges for a0 were estimated to obtain computed D within ±10% of the D used for the generation. Minimum and maximum feature sizes of the profiles were defined and calculated. The feature size range was found to increase with increasing D and σ, in addition to being dependent on the total horizontal length of the profile and the total number of data points in the profile. The suitable range for a0 was found to depend on both D and σ, and then, in turn, on the feature size range, indicating the importance of calculating feature size range for roughness profiles to obtain accurate estimates for the fractal parameters. Procedures are given to estimate the suitable a0 range for a given natural rock joint profile to use with the line scaling method in estimating fractal parameters within ±10% error. Results indicate the importance of removal of global trends of roughness profiles to obtain accurate estimates for the fractal parameters. The parameters C and D are recommended to use with the line scaling method in quantifying stationary roughness. In addition, one or more parameters should be used to quantify the nonstationary part of roughness, if it exists. The estimated C was found to depend on both D and σ and seems to have potential to capture the scale effect of roughness profiles.
 Shirono, T., & Kulatilake, P. H. (1997). Accuracy of the spectral method in estimating fractal/spectral parameters for selfaffine roughness profiles. International Journal of Rock Mechanics and Mining Sciences, 34(5), 789804.More infoAbstract: Accurate quantification of roughness is important in modeling strength, deformability and fluid flow behaviors of rock joints. Selfaffine fractals seem to have the potential to represent rock joint roughness profiles. Both stationary and nonstationary fractional Brownian profiles (selfaffine profiles) with known values of fractal dimension, D, and input standard deviation, a, were generated at different generation levels. A few smoothing techniques were used with the spectral method to calculate D, and two other spectral parameters K, (a proportionality constant; see the text for the details) and CD (the crossover dimension of the profile) for the fractional Brownian profiles. The effects of smoothing, generation level of the profile, seed value used in the generation, nonstationarity of the profile and a on the accuracy of the calculated D were examined using the spectral method. The following conclusions were obtained: (a) To obtain accurate estimates for D, Ks and CD, it seems necessary to have at least W data points per unit length for a profile having a total length of 100 units (this is equivalent to a generation level of 10). (b) For accurate estimation of D, K, and CD, the nonstationarity of profiles should be removed, if it exists, (c) The parameter combinations D and K, (which has the potential to capture scale effects), and D and CD are recommended for quantification of stationary roughness; in addition, extra parameters are required to quantify the nonstationarity. (d) Both the Parzen and Banning smoothing techniques seem suitable to use with the spectral technique to obtain accurate estimates for D, Ks and CD. (e) To obtain accurate estimates for D, Ks, and CD, it is necessary to use a suitable bandwidth for the Parzen window and a suitable number of interations for the Hanning window; this paper provides guidelines to choose these suitable values, (f) Seed value has negligible effect on the accuracy of estimated D, Ks and CD. ©1997 Elsevier Science Ltd.
 Kulatilake, P. H., Chen, J., Teng, J., Shufang, X., & Pan, G. (1996). Discontinuity geometry characterization in a tunnel close to the proposed permanent shiplock area of the Three Gorges dam site in China. International Journal of Rock Mechanics and Mining Sciences and Geomechanics, 33(3), 255277.More infoAbstract: Discontinuity data from a 400m tunnel located close to the proposed shiplock area were investigated to characterize the discontinuity geometry of the rock mass around the tunnel. Traces of 39 major discontinuities (faults and dykes) were found on the tunnel exposures. Over 2000 minor discontinuity (joint) trace data available showed that the rock mass can be separated into about five statistically homogeneous regions. Three to four joint sets were found to exist in each of these regions. Available theoretical probability distributions were found to be insufficient to represent the statistical distribution of orientation of joint clusters. For about 94% of the joint sets, the gamma distribution was found to be the best probability distribution for representing a joint size distribution. Exponential and gamma distributions were the best probability distributions for representing joint spacing distributions for the joint data studied here. Censored joint trace data showed a significant effect of the trace length biases on the estimated mean trace length. This indicated that the 2m wide exposures which were used for collecting the joint trace data may not be sufficient to produce reliable estimates for joint size parameters for the shiplock area. Different estimates were obtained based on wall data and roof data for the joint size and 3D intensity parameters. The following may have contributed to this result: (a) the effect of different finite sample sizes used; (b) the effect of sampling biases; and (c) the effect of modelling uncertainties on the estimated parameters. To obtain better estimates, it is suggested to collect joint data for several exposures which are at least 45m wide and have different orientations. When such data become available for the Three Gorges region, it is suggested to perform validation studies to check the applicability of the models, in addition to performing the joint geometry modelling. Copyright © 1996 Elsevier Science Ltd. All rights reserved.
 Kulatilake, P. H., Shou, G., & Huang, T. H. (1995). Spectralbased peakshearstrength criterion for rock joints. Journal of Geotechnical Engineering  ASCE, 121(11), 789796.More infoAbstract: A new strength criterion, which includes three roughness parameters (D,I, and a spectral coefficient, Ks), is suggested for modeling the anisotropic peak shear strength of rock joints at low normal effective stresses. A validation study shows that the new equation has good capability to predict the anisotropic peak shear strength of joints. In practice, to allow for modeling uncertainties, the new equation should be used along with a factor of safety of about 1.5. Comparisons are also shown between the predicted shear strengths based on Barton's equation and the measured peak shear strengths. from Authors
 Kulatilake, P. H., Shou, G., Huang, T. H., & Morgan, R. M. (1995). New peak shear strength criteria for anisotropic rock joints. International Journal of Rock Mechanics and Mining Sciences and, 32(7), 673697.More infoAbstract: In general, roughness profiles of rock joints consist of nonstationary and stationary components. At the simplest level, only one parameter is sufficient to quantify nonstationary joint roughness. The average inclination angle I, along with the direction considered for the joint surface, is suggested to capture the nonstationary roughness. Most of the natural rock joint surface profiles do not belong to the self similar fractal category. However they may be modelled by selfaffine fractals. Using a new term called specific length, it is shown that even though the fractal dimension D is a useful parameter, it alone is insufficient to quantify the stationary roughness of nonself similar profiles. Also, it is shown why contradictory results for the estimation of D of nonself similar profiles appear in the literature. To estimate D accurately for nonself similar profiles, it seems necessary to use scales of measurement less thanthe crossover length of the profile. Because the crossover dimension of joint roughness profiles can be extremely small, in practice it may be quite difficult to measure roughness at scales of less than the crossover dimension and thus to estimate D accurately. To overcome the aforementioned problems, it is suggested to combine D with a parameter which is negatively correlated to D and also has the potential to compensate for the errors caused by an inaccurate D, and to use the combined parameter to quantify stationary roughness in practice. Four new strength criteria which take the following general form are suggested for modelling the anisotropic peak shear strength of rock joints at low normal effective stresses (00.4 times unconfined compressive strength): τ = σ tanφ + a(SRP)clog10 σJ σd + I where σ, τ, σJ, φ and SRP denote, respectively, the effective normal stress on the joint, peak shear strength, joint compressive strength, basic friction angle, and the stationary roughness parameter. The following four options are suggested to represent the term a(SRP)c: az2′c, aKdbDc, aKsbDc or aKνbDc. Joint roughness data should be used to estimate the roughness parameters I, z2′, kd, Ks, Kν and D in different directions on the joint surface. Parameter D reflects the rate of change in length in response to a change in the scale of measurement r. Because z2′, Kd, Ks and Kν are scaledependent parameters, they can be used to model the scale effect. The coefficients a, b, c and d in the strength criteria should be determined by performing regression analysis on experimental shear strength data. In practice, to allow for modelling uncertainties, the new equations should be used with a factor of safety of about 1.5. © 1995.
 Panda, B., & Kulatilake, P. H. (1995). Study of effect of joint geometry parameters on the permeability of jointed rock. Proceedings of Engineering Mechanics, 1, 337340.More infoAbstract: The equivalent continuum approximation for jointed rock masses was found to be strongly dependent on the relative orientation of joint sets. A strong correlation was found between the directional coefficient of permeability and the fracture tensor component for the connected joint elements in a direction perpendicular to the direction of permeability.
 Kulatilake, P. H. (1994). Scale effects on rock mass deformability. Geomechanics 93. Proc. conference, Ostrava, 1993, 151158.More infoAbstract: Due to the presence of joints, jointed rock masses show anisotropic and scale (size) dependent mechanical properties. A brief description of a recently suggested technique which has the capability of estimating the anisotropic, scaledependent deformability properties of jointed rock is given. Based on this technique, the following have been obtained: 1) 3D plots to show the variation of deformability parameters of jointed rock with joint geometry parameters, 2) a relation between the deformability properties of jointed rock and fracture tensor parameters, 3) an incrementally linear elastic, orthotropic constitutive model to represent the prefailure mechanical behaviour of jointed rock, and 4) some insight to estimate the representative elementary volumes (REVs) and REV property values with respect to the deformability properties of jointed rock. from Author
 Kulatilake, P. H., & Swoboda, G. (1994). Geomechanical modelling of jointed rock. Felsbau, 12(6), 387394.More infoAbstract: A good understanding of the mechanical properties of discontinuous (jointed) rock masses is vital to arrive at safe and economical designs for geotechnical structures in and on rock masses. Geotechnical structures in jointed rock are encountered in civil, mining and petroleum engineering disciplines. The presence of discontinuities strongly affects the mechanical behaviour of rock masses. The following types of problems need to be addressed: a) effect of major discontinuities on deformation and stability of rock masses, b) estimation of equivalent mechanical behaviour of rock masses which contain intact rock and minor discontinuities, and c) procedures to perform stress, deformation and stability analysis of large rock masses which contain intact rock and both major and minor discontinuities. This paper discusses briefly the procedures available, as well as some new techniques, for tackling these problems. from Authors
 Kulatilake, P. H., Ucpirti, H., & Stephansson, O. (1994). Effects of finitesize joints on the deformability of jointed rock at the twodimensional level. Canadian Geotechnical Journal, 31(3), 364374.More infoAbstract: A numerical decomposition technique, which has resulted from a linking between jointgeometry modeling and generation schemes, and a distinct element code (UDEC), is used to study the effect of jointgeometry parameters of finitesize joints on the deformability properties of jointed rock at the twodimensional (2D) level. The influence of jointgeometry parameters such as joint density, ratio of joint size to block size, and joint orientation on the deformability of jointed rock is shown. Relations are established between deformability properties of jointed rock and fracturetensor parameters. An incrementally linear elastic, anisotropic constitutive model is developed to represent the prefailure mechanical behaviour of jointed rock at the 2D level. This constitutive model has captured the anisotropic, scaledependent behaviour of jointed rock. In this model, the effect of the jointgeometry network in the rock mass is incorporated in terms of fracturetensor components. Some insight is given related to estimation of representative elementary volumes for deformability properties of jointed rock.
 Kulatilake, P. H., Wang, S., & Stephansson, O. (1993). Effect of finite size joints on the deformability of jointed rock in three dimensions. International Journal of Rock Mechanics and Mining Sciences and, 30(5), 479501.More infoAbstract: A numerical decomposition technique, which has emerged from a linking between joint geometry modelling and generation schemes and a distinct element code (3DEC), is used to evaluate the effects of joint geometry parameters of finite size joints on the deformability properties of jointed rock at the threedimensional (3D) level. Variation of deformability parameters of jointed rock with joint geometry parameters such as joint density, joint size/block size and joint orientation, are shown through 3D plots. Relations are developed between deformability properties of jointed rock and fracture tensor parameters. An incrementally linear elastic, orthotropic constitutive model is suggested to represent the prefailure mechanical behaviour of jointed rock. This constitutive model has captured the anisotropic, scaledependent behaviour of jointed rock. In this model, the effect of the joint geometry network in the rock mass is incorporated in terms of fracture tensor components. Some insight is given related to estimation of Representative Elementary Volumes (REVs) and REV property values with respect to deformability properties of jointed rock. © 1993.
 Kulatilake, P. H., Wathugala, D. N., & Stephansson, O. (1993). Joint network modelling with a validation exercise in Stripa mine, Sweden. International Journal of Rock Mechanics and Mining Sciences and, 30(5), 503526.More infoAbstract: In this paper, eight joint geometry modelling schemes are suggested and applied to a set of Stripa mine joint data to build 3D joint networks to a granitic rock mass. These modelling schemes include investigations for statistical homogeneity of the rock mass, corrections for sampling biases, and applications of stereological principles, to estimate 3D joint geometry parameters from 1D or 2D joint geometry parameter values. Results show the possibility of obtaining different estimates for both joint size and joint intensity parameters through these different schemes. This indicates the importance of performing validation studies for developed joint geometry modelling schemes. Validation procedures developed and performed indicated the need to try out different schemes in modelling joint geometry parameters in order to establish realistic 3D joint geometry modelling schemes which provide good agreement with field data during verification. It is important to realize that different types of joint geometry modelling schemes are needed to model joint networks in different types of rock formations. © 1993.
 Kulatilake, P. H., Wathugala, D. N., & Stephansson, O. (1993). Stochastic three dimensional joint size, intensity and system modelling and a validation to an area in Stripa Mine, Sweden. Soils and Foundations, 33(1), 5570.More infoAbstract: In this paper, eight joint geometry modelling schemes are suggested and applied to a set of Stripa mine joint data to build threedimensional joint networks to a granitic rock mass. These modelling schemes include investigations for statistical homogeneity of the rock mass, corrections for sampling biases, and applications of stereological principles, to estimate 3D joint geometry parameters from 1D or 2D joint geometry parameter values. Results show the possibility of obtaining different estimates for both joint size and joint intensity parameters through these different schemes. This indicates the importance of performing validation studies for developed joint geometry modelling schemes. Validations performed indicated the need to try out different schemes in modelling joint geometry parameters in order to come up with realistic 3D joint geometry modelling schemes which provide good agreement with field data during verification. It is important to realize that different types of joint geometry modelling schemes are needed to model joint networks in different types of rock formations.
 Wang, S., & Kulatilake, P. H. (1993). Linking between joint geometry models and a distinct element method in three dimensions to perform stress analyses in rock masses containing finite size joints. Soils and Foundations, 33(4), 8898.More infoAbstract: To use the distinct element method to perform stress analysis in three dimensions (3D), the domain considered has to be discretized into distinct polyhedral blocks by the joints present in rock masses. Unfortunately, most of the joint networks present in rock masses do not produce distinct polyhedral blocks. Therefore, to use the distinct element method to perform stress analysis in actual jointed rock masses, it is necessary to find a procedure to discretize the rock mass into polyhedral blocks starting from a 3D joint network which does not discretize a rock mass into polyhedral blocks. Also, this procedure should not alter the mass deformation and strength properties of the jointed rock. This paper provides needed procedures to link 3D stochastic or deterministic joint geometry modeling schemes to a 3D distinct element method to perform stress analysis in jointed rock masses which contain finite size joints. The procedure is illustrated through an example.
 Kulatilake, P. H., Ucpirti, H., Wang, S., Radberg, G., & Stephansson, O. (1992). Use of the distinct element method to perform stress analysis in rock with nonpersistent joints and to study the effect of joint geometry parameters on the strength and deformability of rock masses. Rock Mechanics and Rock Engineering, 25(4), 253274.More infoAbstract: To use the distinct element method, it is necessary to discretize the problem domain into polygons in two dimensions (2 D) or into polyhedra in three dimensions (3 D). To perform distinct element stress analysis in a rock mass which contains nonpersistent finite size joints, it is necessary to generate some type of fictitious joints so that when they are combined with the nonpersistent joints, they discretize the problem domain into polygons in 2 D or into polyhedra in 3 D. The question arises as to which deformation and strength parameter values should be assigned to these fictitious joints so that they behave as intact rock. In this paper, linear elastic, perfectlyplastic constitutive models with the MohrCoulomb failure criterion, including a tension cutoff, were used to represent the mechanical behaviour of both intact rock and fictitious joints. It was found that, by choosing the parameter values for the constitutive models as given below, it is possible to make the fictitious joints behave as intact rock, in a global sense. a) For both the intact rock and the fictitious joints, the same strength parameter values should be used. b) A joint shear stiffness (JKS) value for fictitious joints should be chosen to produce a shear modulus/JKS ratio (G/JKS) between 0.008 and 0.012 m. c) A joint normal stiffness/JKS ratio (JKN/JKS) between 2 and 3 should be chosen. The most appropriate value to choose in this range may be the Young's modulus/G value (E/G) for the particular rock. Some examples are given in the paper to illustrate how to use the distinct element method to perform stress analysis of rock blocks which contain nonpersistent joints and to study the effect of joint geometry parameters on strength and deformability of rock masses. © 1992 SpringerVerlag.
 Kulatilake, P. H., Wang, S., & Uepirti, H. (1992). Joint network modeling and scale effects in rock. Proceedings of Engineering Mechanics, 441444.More infoAbstract: Joint networks play a vital role in the mechanical and hydraulic behavior of rock masses. A realistic representation of joint geometry in rock masses and correct understanding of the mechanical properties of rock masses are very important in solving rock engineering problems. The first part of this paper shows how to build threedimensional stochastic joint network schemes and how to perform validation studies for these schemes. An investigation and the results obtained so far on the effects of joint geometry parameters on the mechanical properties of jointed rock masses are presented in the second part of the paper.
 Kulatilake, P. H., & Lacasse, S. (1991). Probabilistic equivalent linear soil spring stiffness analysis for gravity platforms: Conceptual model. Computers and Geotechnics, 12(1), 128.More infoAbstract: Horizontal and rotational soil spring stiffnesses of the foundation are needed (a) to perform the dynamic analysis of gravity platforms and (b) to obtain characteristic wave loads to be used in the geotechnical design analysis of platforms. Dynamic analyses are especially important in the design of deepwater platforms. A probabilistic procedure in the context of the first two moments is suggested to estimate the equivalent linear soil spring stiffnesses of the foundation for gravity platforms on clay loaded under undrained conditions. In calculating the first two moments of the soil spring stiffnesses, consideration is given to the influence of the uncertainties in the loads and soil properties, and the model errors encountered in applying the deterministic procedure. Limitations in the available data, the current analytical procedures, and in the probabilistic analysis are discussed. The companion paper [1] provides an illustrative example on the same topic. © 1991.
 Kulatilake, P. H., & Wathugala, D. (1991). Stochastic threedimensional joint geometry modelling including a verification to an area in Stripa Mine, Sweden. Array, 545555.More infoAbstract: At present, a threedimensional joint geometry modelling scheme which investigates statistical homogeneity, incorporates corrections for sampling biases and applications of stereological principles, and, also includes a formal verification procedure is not available in the literature. This paper shows development of 3D joint geometry models with aforementioned features. Joint data from Stripa mine were used in the investigation. Verifications performed so far has indicated the need to try out different schemes in modelling joint geometry parameters and to consider a number of combinations of such schemes in order to arrive at a realistic 3D joint geometry model which provides a good comparison with field data during verification. Further verification studies are recommended to check the validity of the suggested modelling schemes.
 Kulatilake, P. H., Lacasse, S., & Gabr, M. (1991). Probabilistic equivalent linear soil spring stiffness analysis for gravity platforms: Illustrative example. Computers and Geotechnics, 12(1), 2954.More infoAbstract: The companion paper [1] suggests a procedure to perform equivalent linear soil spring stiffness analysis for gravity platforms founded on clay and loaded under undrained conditions. This paper provides an example to illustrate the suggested procedure. The CONDEEP SP41024 platform concept proposed for the Troll field was used for the example. Plastic Drammen clay was chosen as the foundation soil to a depth of 150 m. Wave loads on the platform were estimated assuming a 70 yr platform life, plus an 18 hr storm buildup, plus a 6 hr storm with a return period of 100 yrs. The constitutive model for the soil was calculated for the 6 hr rare event in the load history, assuming that the rare even occurs at the end of the load history. The coefficients of variation of the equivalent soil spring stiffnesses (horizontal and rotational) were found to be around 0.50. © 1991.
 Kulatilake, P. H., Wu, T. H., & Wathugala, D. N. (1990). Probabilistic modelling of joint orientation. International Journal for Numerical and Analytical Methods in Geomechanics, 14(5), 325350.More infoAbstract: Observed frequencies of joint orientations are subject to error due to sampling bias. This error should be corrected before statistical inference is made on the distribution of orientation. Corrections (weighting functions) are developed for sampling bias in orientation for finite joints of different sizes and shapes intersecting rectangular exposures. Chisquare goodnessoffit procedures available for hemispherical normal and bivariate normal distributions are modified to make them applicable for both raw data and data corrected for sampling bias. The aforementioned corrections and procedures were applied to a joint orientation cluster to study the effect of (a) joint orientation, (b) joint size and (c) joint shape, on the statistical distribution of the orientation. The influences of all these aforementioned factors were found to be significant. However, at present, joint sizes and shapes are not measured in field joint surveys. Therefore, it is suggested to make attempts to obtain joint sizes and shapes in field joint mapping surveys.
 Wathugala, D. N., Kulatilake, P. H., Wathugala, G. W., & Stephansson, O. (1990). A general procedure to correct sampling bias on joint orientation using a vector approach. Computers and Geotechnics, 10(1), 131.More infoAbstract: Observed relative frequencies of joints should be corrected for sampling bias before inferring statistical distributions for orientations. The procedure available for sampling bias correction when finite size joints intersect finite size sampling domains is directly applicable only for vertical sampling planes [13]. Thus, a general procedure applicable for sampling domains of any orientation is developed in the present study. The corrected frequency of a joint is obtained by assigning a weight to each joint through a weighting function which is inversely propotional to the probability of intersection between the joint and the sampling plane. This probability of intersection is determined from a hypothesis indicating "the probability of intersection is proportional to the volume in which the center of the sampling domain should lie in order to intersect the joint". A vector approach to find this volume is described herein, followed by applications of this method to study the influence of the sampling bias correction on orientation frequency and to find statistical distributions for orientation data. © 1990.
 Ghosh, A., & Kulatilake, P. H. (1989). Reply to correction for "a FORTRAN program for generation of multivariate normally distributed random variables". Woronow, A., 1989, Computers & Geosciences, v. 15, no. 6, p. 1033. Computers and Geosciences, 15(6), 10341035.
 Kulatilake, P. H., & Varatharajah, P. (1989). A weighted regression approach to model and estimate spatial variability of soil properties in one dimension. Computer and physical modelling in geotechnical engineering. Proc. international symposium, Bangkok, 1986, 279285.More infoAbstract: A method is developed to estimate spatial variability of soil properties in one dimension. An example is given to illustrate the use of the developed technique. With slight modifications, the method can be extended to estimate spatial variation in three dimensions. from Authors
 Ghosh, A., & Kulatilake, P. H. (1987). A FORTRAN program for generation of multivariate normally distributed random variables. Computers and Geosciences, 13(3), 221233.More infoAbstract: The computer program given in this paper generates a set of values for each of the random variables which are distributed according to a multivariate normal distribution. It is written in FORTRAN 77 and is designed to run on a CYBER 175 computer. In generating a set of values, the program either can use actual data of the variables as input to estimate parameter values of the multivariate normal distribution or the parameter values of the multivariate normal distribution can be used directly as input to the program. The theory and the necessary algorithms for the generation are given in detail. Use of the program is illustrated through an example in soil engineering. MonteCarlo simulation method is used for working out the example. © 1987.
 Kulatilake, H. S., & Miller, K. M. (1987). SCHEME FOR ESTIMATING THE SPATIAL VARIATION OF SOIL PROPERTIES IN THREE DIMENSIONS.. Array, v 2p.More infoAbstract: A procedure is developed to estimate the spatial variability of soil properties in 3D from irregularly spaced data. Two basic components of the soil property data are modeled: (1) the global trend which is nonstationary and is estimated using multiple linear regression analysis for the raw soil property data; and (2) the local trend which contains a stationary portion that is characterized by its autocorrelation structure. Autocorrelation between residual values is determined and the correlation structure is modeled using multiple linear regression analysis. Local trend is estimated using a kriging process that utilizes the autocorrelation structure. The estimation model was applied to a set of Dutchcone penetration test data. Analysis on the residual values showed that autocorrelation was weak, which resulted in a large estimated variance for the estimated mean property value. It was concluded that better soil property estimates could be obtained using smaller sampling intervals. The model can also be used in the design of soil exploration programs.
 Kulatilake, P. H. (1987). Modelling of cyclical stratigraphy using Markov chains. International Journal of Mining and Geological Engineering, 5(2), 121130.More infoAbstract: Stateoftheart on modelling of cyclical stratigraphy using firstorder Markov chains is reviewed. Shortcomings of the presently available procedures are identified. A procedure which eliminates all the identified shortcomings is presented. Required statistical tests to perform this modelling are given in detail. An example is given to illustrate the presented procedure. © 1987 Chapman and Hall Ltd.
 Wu, T. H., Williams, R. L., Lynch, J. E., & Kulatilake, P. H. (1987). STABILITY OF SLOPES IN RED CONEMAUGH SHALE OF OHIO.. Journal of geotechnical engineering, 113(3), 248264.More infoAbstract: Nine slopes excavated in soft Conemaugh shales containing slickensides were investigated. The trace lengths and orientations of the slickensides found in exposures were measured. The shear strengths of the intact shales and slickensides were measured by means of direct shear tests. Stability analyses were performed to calculate the shear strength of the shale in the failed slopes. For the stable slopes, the strength required to attain a safety factor of 1 was calculated. The in situ strength must exceed this strength. The estimated strengths are examined in the light of the average density and the average length of slickenside traces and the stress distribution. Estimates are made of the effect of slickensides and stress distribution on the reduction of shear strength below its peak value.
 Kulatilake, P. H. (1986). Bivariate normal distribution fitting on discontinuity orientation clusters. Mathematical Geology, 18(2), 181195.More infoAbstract: A bivariate normal density function has been used to represent discontinuity orientation cluster distributions. Goodnessoffit tests should be performed in order to make decisions on the representation of discontinuity clusters by theoretical probability distributions. In the literature, graphical procedures are available to fit a bivariate normal distribution to discontinuity clusters. However, these procedures assume no correlation between the two orientation parameters. In this paper (a) a numerical procedure, and (b) a semigraphical procedure are given to perform a κ{script}2 goodnessoffit test for bivariate normal distributions having nonzero correlation coefficient between the two parameters. These procedures were applied to a selected discontinuity cluster. The semigraphical procedure was found to be a timeconsuming process. On the other hand, rapid computation can be done with the computer program developed for the numerical method. Sensitivity of the κ{script}2 test results of the IXJ grid setup was investigated. Mean orientation estimation for the cluster based on the equal area polar projection was compared with the estimation based on the moment estimate method. For the cluster analyzed, estimations of bivariate normal parameters {Mathematical expression} and {Mathematical expression}, based on the equalarea polar projection, and values based on the moment estimation method were found to be different up to about 6.7% of the values based on the moment estimation method. © 1986 Plenum Publishing Corporation.
 Kulatilake, P. H., & Fuenkajorn, K. (1986). CALCULATION OF MINIMUM ROCK BOLT FORCE AND MINIMUM STATIC ACCELERATION IN TETRAHEDRAL WEDGE STABILITY.. Array, 197203.More infoAbstract: Hoek, Bray and Boyd (1973) have given an analytical procedure to calculate the factor of safety for the case of potential sliding along the line of intersection of the two supporting planes for a threedimensional wedge having two free surfaces. In this paper, this analytical procedure is extended to obtain expressions to calculate minimum rock bolt force and minimum static acceleration to achieve required factor of safety for the same type of wedge stability. The effects of an obliquely inclined upper slope face, a tension crack behind the slope face, water pressure within the tension crack and on the two supporting planes and any other known external forces can be included in the analysis. Resistance to sliding results from both friction and cohesion on the two supporting surfaces. Examples are given to illustrate the use of derived expressions.
 Kulatilake, P. H., & Fuenkajorn, K. (1986). Calculation of minimum rock bolt force and minimum static acceleration in tetrahedral wedge stability.. Geotechnical stability in surface mining. Proc. symposium, Calgary, 1986, 197203.More infoAbstract: An analytical procedure is extended to obtain expressions to calculate minimum rock bolt force and minimum static acceleration to achieve required factor of safety for the same type of wedge stability. The effects of an obliquely inclined upper slope face, a tension crack behind the slope face, water pressure within the tension crack and on the two supporting planes and any other known external forces can be included in the analysis. from Authors
 Kulatilake, P. H., & Ouyang, S. (1986). PROBABILISTIC MODELING OF SHEAR STRENGTH OF ROCK DISCONTINUITIES USING DIRECT SHEAR TEST DATA.. Proceedings  Symposium on Rock Mechanics, 112120.More infoAbstract: Available literature on the topic is reviewed. The MohrCoulomb model and a power law model are used as strength criteria to model shear strength of rock discontinuities. Use of various regression techniques for point estimation of shear strength parameters from direct shear data are discussed. Bivariate beta and bivariate normal distributions have potential to represent the statistical distribution of strength parameters.
 Kulatilake, P. H., & Wu, T. H. (1986). RELATION BETWEEN DISCONTINUITY SIZE AND TRACE LENGTH.. Proceedings  Symposium on Rock Mechanics, 130133.More infoAbstract: The concepts of conditional probability and geometrical probability are used to establish a relation between the probability density of the trace length and the probability density of the discontinuity diameter. The discontinuity is represented as a thin circular disc. An example is used to illustrate the possibly large difference between the mean discontinuity diameter and the mean trace length.
 Kulatilake, P. H. (1985). ESTIMATING ELASTIC CONSTANTS AND STRENGTH OF DISCONTINUOUS ROCK.. Journal of Geotechnical Engineering, 111(7), 847864.More infoAbstract: In this paper, a numerical method that combines analytical decomposition modeling with statistical simulation is presented. The method consists of: Probabilistic modeling of discontinuity geometry; generation of discontinuities in rock blocks by MonteCarlo simulation; and finite element analysis of simulated rock blocks. The method was used to study the mass properties of a shale that contains slickensides. For ten simulations, the coefficient of variations of all three mass parameters were found to be less than 0. 21.
 Kulatilake, P. H., & Wu, T. H. (1984). Estimation of mean trace length of discontinuities. Rock Mechanics and Rock Engineering, 17(4), 215232.More infoAbstract: Trace lengths of discontinuities observed on finite exposures are biased due to sampling errors. These errors should be corrected in estimating mean trace length. A technique, which takes into account the sampling errors, is proposed for estimating the mean trace length on infinite, vertical sections from the observations made on finite, rectangular, vertical exposures. The method is applicable to discontinuities whose orientation is described by a probability distribution function. The method requires that the numbers of discontinuities with both ends observed, one end observed, and both ends censored be known. The lengths of observed traces and the density function of trace length are not required. The derivation assumes that midpoints of traces are uniformly distributed in the vertical plane. Also independence between trace length and orientation is assumed. Data on a Pennsylvania shale in Ohio, U. S. A., were used as an example. © 1984 SpringerVerlag.
 Kulatilake, P. H., & Wu, T. H. (1984). SAMPLING BIAS ON ORIENTATION OF DISCONTINUITIES.. Rock mechanics and rock engineering, 17(4), 243253.More infoAbstract: In discontinuity surveys, biases are frequently introduced in sampling for geometrical properties. Terzaghi developed a correction for orientation data for planar discontinuities of infinite size. For discontinuities of finite size intersecting finite size exposures, the shape and size of both the discontinuity and the exposure also influence the probability of intersection. In discontinuity modeling, joints have been represented as thin circular discs and as Poisson planes. Robertson reported that strike length and dip length of discontinuities were approximately equal. In this paper the discontinuities are considered to be thin circular discs with a finite diameter. The probability of the disc of a given diameter intersecting a vertical plane is derived and used to correct observed orientation data.
 Kulatilake, P. H., & Wu, T. H. (1984). The density of discontinuity traces in sampling windows. International Journal of Rock Mechanics and Mining Sciences and, 21(6), 345347.