A New Approach to Rapid 3D Mapping of Rock Mass Structure
The 12th International Conference of
International Association for Computer Methods and Advances in Geomechanics (IACMAG)
1-6 October, 2008
A New Approach to Rapid 3D Mapping of Rock Mass Structure
Seydisehir Voc. Sch. of High. Educ., Selçuk University, Konya, Turkey
Department of Computer Engineering, Selçuk University, Konya, Turkey
Key Words: Discontinuity, rock mass, mathematical transformation, isometric perspective, 3D mapping.
ABSTRACT: The prediction of rock mass behavior is an important task in many engineering projects, as the
behavior of rock masses can be controlled by the presence of discontinuities. Being able to map the structure
of a rock mass is crucial to understanding its potential behavior. This understanding can positively impact the
safety and efficiency of an engineering project. In this research, rock masses were mapped and analyzed using
linear mathematical transformations and isometric perspective methods to achieve meaningful three-dimensional
(3D) results. The rock mass fracture representation is based on explicit mapping of rock faces. The developed
model can improve safety and productivity through its application in the determination and analysis of rock
mass structure for geological and geotechnical assessment. Based on the methods explained here, a software
system was developed for analyzing the geometric characteristics of discontinuities in a rock mass. In this model,
discontinuities in a rock mass can be visualized both individually and in combination, and cross-sections can be
generated at any desired location. In addition, intersection lines between discontinuities can be generated as dip
direction vectors. The natural structure attained by using this developed model agrees well with field
Visualization is the task of generating and understanding images that contain important features. Visual sight
constitutes 70% of sense of object perception (Ming and Peter, 2006), which is valid in engineering applications,
where the structural reconstruction of an actual object is a required step for accurate visualization. To
comprehend, render, and reveal complex structural objects such as in open pit mines, appropriate geometric
models must be designed. Two-dimensional (2D) models, including geological maps, cross-sections, sketches of
strain and stress patterns, and stereographic nets, are used in mining operations. Normally, the set of
observations and measurements supporting these models is small in relation to the complexity of the real objects
from which they are derived. Therefore, these models need additional guidance and explanations. However,
geological modeling techniques have evolved from highly conceptual methods to practical computing methods,
which include three-dimensional (3D) approaches (Zhong et al., 2006). 3D models are accepted as promising
tools for a better comprehensive understanding of engineering object characteristics. Accurate, detailed 3D
representations of real objects are needed for successful engineering design. In addition, suitable
representations can prove to be useful as strong communication tools between technical experts and nonspecialists.
In designing large and risky engineering structures, it is advised to first conduct research on models to determine
alternative solutions before deciding on the best alternative, and in order to reveal the optimum economical
boundaries and the feasibility of the project. For these aims, one of the modeling methods used for technological
purposes is geometric modeling, in which the similarities between the model and a scaled prototype of the model
are considered. In other words, there is a fixed ratio between the corresponding points for any specific feature in a
coordinate system. In rock masses, these models are made from discontinuity intersection lines as 2D and 3D
networks. In many geomechanical models, the most commonly used geometrical features are spacing,
orientation, and trace length of discontinuities, which have been defined in detail previously (ISRM, 1978).
The discrete element method (DEM) is a widely used modeling technique in rock engineering applications. The
foundation of the DEM was first developed by Cundall (1971). Several researchers (Hocking, 1977; Williams et
al., 1985; Pande et al., 1990; Shi, 1996; Jing, 1998) then introduced additional aspects of the DEM in different
related engineering problems. The key concepts of the DEM are that the domain of interest is treated as an
assemblage of rigid or deformable blocks/particles, and that the contacts between them need to be identified and
continuously updated during the entire deformation/motion process and be represented by appropriate
constitutive models. The discrete fracture network (DFN) approach has been also widely used to investigate the
behavior of rock masses, especially the hydraulic behavior of fractured rock masses. This technique was
developed in the early 1980s for both 2D and 3D problems (Long et al., 1982; Andersson, 1984; Endo et al.,
1984; Smith and Schwartz, 1984; Elsworth, 1986; Dershowitz and Einstein, 1987). The DFN method considers
fluid flow and transport processes in fractured rock masses through a system of connected fractures. There are a
number of software programs derived from DEM and DFN which are available. The most representative explicit
DEM methods are the computer codes UDEC and 3DEC (by Itasca Consulting Group Inc.) for 2D and 3D
problems in rock mechanics. The most widely known DFN codes are NAPSAC (by AEA Decommissioning &
Radwaste, Harwell, UK) and FRACMAN/MAFIC (by Golder Associates, Inc.) for 2D and 3D problems in rock
mechanics. These software programs have been successfully used to simulate the behavior of rock fractured
systems. However, their application to engineering problems is still limited to simplified representations of rock
masses, and the inclusion of discontinuities is still a challenging task.
In this 3D modeling research, we aimed to provide accurate and practical information and details for engineering
applications such as rock stability, pre-design of mines, blasting, and connectivity of fluids. In the model,
geometric features of discontinuities measured with scan-line surveys on the rock mass outcrop, such as dip, dip
direction, and spacing, were used as basic geometric features of the rock mass. The main instruments for taking
data were a compass-clinometer and a tape meter. Additionally, mathematical transformations, perspective rules
and developed schemes, and mathematical approaches and equations were used for the developed model.
Moreover, linear relations were used and an isometric presentation was preferred for all illustrations.
2 The basics of the proposed 3D modelling method
The modeling sequence followed includes: a) identification of spatial positions of the discontinuity planes and the
outcrop as four types; b) mathematical definition of the discontinuity planes; c) isometric transformation and
projection of the rock structure, including discontinuity planes and intersection lines; d) generation of the
intersection lines as dip vectors; and e) developing cross-sections from the obtained 3D model.
2.1 Classification of discontinuities positions
In this study, the first step is the classification of the discontinuity orientations. The classification term is a term
used to define the reorganizing of the classical discontinuity orientations into classes. Discontinuity and rock mass
outcrop orientations may be different in nature relative to north. Because of this, two strategies can be developed
in the model. The first strategy uses four possible alternative discontinuity orientations in a coordinate system
(Figure 1(a)). In the boundary values of the discontinuity orientation, the apparent dip, dip direction, and outcrop
(YZ plane) orientations are considered. The angles used in the calculation are shown in Figure 1(b). This figure
was drawn for the first discontinuity type shown in Figure 1 (a). The second strategy uses 64 possible
discontinuity orientations according to four different possible discontinuities. (Figure 1(c)). In this figure, the
possible spatial orientations of the discontinuities are given on the XY plane.
Figure 1. Classification of the possible orientations of any discontinuity; (a) four possible general discontinuity
orientations, (b) explanation of the angular conditions and values and numeric results for the 1/1 (column/row)
type discontinuity (α: dip direction, γ: outcrop orientation, the search angle is shown with a dashed line), (c)
possible orientations for these four discontinuities relative to the north (on the XY plane).
The process of classifying discontinuities is the first step of our proposed model. Discontinuity planes and rock
mass outcrops can be naturally in different orientations relative to north. Because of this, the two following
strategies, which depend on each other, were developed. These strategies are related to 3D discontinuity planes
in the rock mass and 2D discontinuity traces on the visible surfaces of the rock mass, respectively. In the first
strategy, the four possible alternative discontinuity orientations in a constant coordinate system and the north
(Figure 1(a)) were constructed. In these illustrations, discontinuities perpendicular and parallel to the XZ, XY, and
YZ surfaces are neglected. The dip angle (β) (on the YZ surface), dip direction angle (α), and the outcrop
orientation angle (γ; on the XY surface) were considered as the boundary values of the discontinuity orientation.
The dip direction, outcrop orientation (in this paper, we used spatial positions of the outcrop relative to north as
the outcrop orientation), and the unknown angle (dashed line) which is used in the calculations are also shown in
Figure 1(b). Here, the outcrop surface was assumed to be vertical, and this surface was taken as the YZ surface.
This figure was generated using the first discontinuity type as an example (details of this type are also shown in
Figure 1). In the second strategy, 16 configurations for each discontinuity type (Figure 1 (c)) were developed
relative to north. Thus, a total of 64 configuration types were obtained for the four discontinuity types. Here, the
XY plane was divided into eight equal areas to assess the possible conditions. In this paper, the angular condition
and angular value definitions are proposed as a new concept.
In Figure 1(c), outcrop orientations were classified relative to north in every row as 16 configurations. In addition,
discontinuity types are classified according to the constant outcrop position in each of the 16 columns for the four
configurations. All these illustrations are on the XY surface. Derived angular conditions and values for the 1/1, 1/2,
1/3 and 1/4 (column/row) discontinuity types (Figure 1 (a and c)) are presented in Table 1. It is possible to write
three or four angular conditions and three angular values for each configuration. Consequently, a total of 192
angular values were obtained. In our method, the discontinuity type (Figure 1(a); type 1, 2, 3, or 4) is selected
first. The final angular values are then obtained by utilizing the angular conditions. For each angular condition
configuration, there is an angular value. An angular value consists of three sub-elements which are on the x, y,
and z axes. Discontinuity type 1/1 (the first type) in Figure 1(b) can be considered as an example of the geometric
approach. Only the β (dip) angle can be used on the YZ plane for the drawing of the discontinuity trace. During
this process, the three angles considered are 90°, β, 270°+β clockwise from the discontinuity trace to the x-, y-,
and z-axes, respectively. On the XY plane, the angular conditions formed should be taken into account. Here, the
x-axis was the North direction, and the required angle represented by a dashed line can be calculated by using
the outcrop orientation angle (γ) and the dip direction angle (α) which was recorded from the scan line survey.
The required angle is taken as clockwise between the discontinuity trace and the x-axis.
Table 1. The angular conditions and values for types 1/1, 1/2, 1/3 and 1/4 (column/row).
2.2 Mathematical definition of the discontinuity planes
After assigning the classes for the discontinuities, the boundary points of the discontinuity will determine the
modeling of the discontinuity in the rock block. Afterwards, the discontinuity plane is generated from these points.
The method used for this is discussed in this section. Let the angles of two straight lines be αi, βi, and γi related to
the x-, y-, and z-axes, respectively. The parametric equations used for calculating the intersection point between
these two straight lines are as follows:
( 1 i1 ) 2 ( i2 ) 3 x + λ * Cos α = x + λ * Cos α = x (1)
( ) ( ) 1 i1 2 i2 3 y + λ * Cos β = y + λ * Cos β = y (2)
( ) ( ) 1 i1 2 i2 3 z + λ * Cos γ = z + λ * Cos γ = z (3)
where the left parts of the equations are the parametric equality of a straight line and the remaining parts are the
parametric equality of another straight line and the coordinates of the intersecting point between two straight lines,
respectively. For every adjacent discontinuity trace as a straight line in the equation system, a λ value and a t
value should be calculated. In these equations, the first equation (Eq. 1) does not fit the solution; the software
reaches the conclusion according to Eqs. 2 and 3. In Eq. 2, one of the unknowns is given in the form of the
others, and the solution is obtained using Eq. 3. The λ value and the t value are calculated using Eqs. 4 and 5.
( ) ( )
( ) ( )
( ) ( )⎥⎦
Cos β * Cos γ
* Cos γ
z z y y
( ) * t
2 1 ⎟
λ = (5)
Using the obtained λ and t values, the intersection point values of the p (x3, y3, z3) can be found and recorded in a
file. In this study, the visible discontinuity traces on the outcrop and prism edges were assumed to be straight
lines. The discontinuity simulation process is based on constructing the edges which formed the discontinuity
plane in a prism. This calculation process is performed on the visible surfaces of the prism. Thus, all intersection
points between edges belonging to prism and the prism edge–discontinuity trace can be obtained using the
aforementioned equations. The discontinuity traces can be drawn on visible surface of the prism as straight lines
using the intersection points which are calculated on the prism edges. These lines are projected to the hidden
edges of the prism. The obtained four straight lines which bounded the discontinuity plane were used as a tool in
the simulation of the discontinuity plane. In addition, the intersection lines between the discontinuity planes can
also be derived using similar processes.
In this paper, an isometric transformation process was applied for all intersection points belonging to the
discontinuity trace–discontinuity trace and the discontinuity trace–prism edge. These points were used in the
software as the construction elements, as explained below.
3 The developed software
In this study, a program was developed to simulate the rock structure. The Visual Basic 6.0 programming
language was used in the software, and the inputs of the package program are the data obtained from the outcrop
of the rock mass by the scan-line method. To make a full representation of each discontinuity, the dip and dip
direction are recorded. An identity number (ID) was given to each discontinuity. Whereas dip and dip direction
data are variable, there are also constant data used to simulate the rock mass, including the dimensions of the
prism in which the rock mass is simulated, the scan-line height, and outcrop orientation (γ) angles between the
rock mass and north. An input table is formed by adding these values to the discontinuity data. The input table is
called the information system table (IST). Each row of the IST represents information related to a discontinuity,
and each column represents the qualities of the rock mass. There are five columns in the IST. These columns are
discontinuity ID, dip direction (α), dip (β), outcrop orientation (γ), and cumulative spacing. An example of the IST
which was used in the field experiment is given in Table 2.
The calculation which enables the simulation of the rock mass can then be initiated. The calculations perform the
following steps in turn. Let us consider discontinuity is ID-1. With the dip and dip direction information, the class of
ID-1 is determined (discontinuity type one, two, three, or four). By determining the class of ID-1, its line on the YZ
plane |AB| is calculated. The points which pass through the YZ plane of the line (A and B) must be determined.
The approach, which gives the intersection of two discontinuities and which is newly suggested by us, is used
here. The first discontinuity is ID-1, whereas the second discontinuity is the line RM1RM2. The intersection of the
two discontinuities A is calculated using Equations 1, 2, and 3. The second discontinuity is called the RM3RM4
line. Point B is obtained using the suggested approach. Here, point RM1 is taken as (0, 0, 0), point RM2 is taken
as (0, Ybd, 0), point RM3 is taken as (0, 0, Zbd), and point RM4 is taken as (0, Ybd, Zbd). The same processes
are used in the calculation of points C and D. Consequently, discontinuity plane ID-1 is obtained (Figure 2). Points
A, B, C, and D are recorded for later use. After the discontinuity plane is obtained, it is tested to determine
whether it crosses over the prism edges. If there is a crossover, the discontinuity plane is trimmed with the prism
boundaries. The software can obtain all discontinuity planes in this way.
Figure 2. The simulation of discontinuity plane ID–1 in the selected prism.
The intersection of the discontinuity planes is a problem that is solved in the software. All intersections can be
detected. If there are detected intersections, the intersection lines are calculated. These detection processes are
dealt after all discontinuity planes are obtained. For instance, discontinuity ID-1 is compared with discontinuities
ID-2, ID-3,…, ID-N on the YZ plane. The plane coordinates of two different discontinuities are found on both the
YZ plane and the XY planes. For both discontinuities, the traditional mathematical approach is used to detect
whether the plane lines are parallel to the YZ plane or the XY plane. If there is no parallelism, there is an
intersection point. With the help of Equations 1, 2, 3, 4, and 5, the intersection point K, which is between the two
discontinuities on the YZ plane, is calculated. Intersection point K can be on the YZ plane of the rock mass, or it
could possibly also be outside the boundaries of the mass. For this reason, the control of point K with the mass
boundaries is also performed. Similarly, point L, which is on the XY plane, is also calculated. These processes are
performed on all planes of the rock mass except for the XZ planes. All obtained intersection points are recorded.
There is a possibility that a discontinuity plane can intersect more than one discontinuity plane. In this case, the
discontinuity is compared with all the remaining planes. These calculations increase the flexibility of the image in
the software, and each discontinuity can be examined individually.
The general flow diagram and a typical interface of the developed software are shown in Figures 3 and 4. As seen
in Figure 3, the building of the information system table includes recording data from the first phase. For the
calculation of a discontinuity plane, four intersection points are computed on the prism edges. Then, the
discontinuity plane is generated from these points. Moreover, the software can find the intersection point for any
two or more discontinuity planes and the corresponding intersection lines. Thus, the images of all discontinuity
planes are simulated.
Table 2. An example of the IST that was used in the field experiment.
Discontinuitiy Dip Dip Outcrop Cumulative Points
number Direction Orientation Spacing (m)
(α) (β) (γ) (m) X Y Z X Y Z
1 160 113 23 2,15 0 2,15 1,5 8 35 4
2 246 55 23 2,6 0 2,6 1,5 8 35 4
3 249 115 23 9,15 0 9,15 1,5 8 35 4
4 267 110 23 10,05 0 10,1 1,5 8 35 4
5 193 52 23 12,05 0 12,1 1,5 8 35 4
6 250 127 23 17,1 0 17,1 1,5 8 35 4
7 252 132 23 20,85 0 20,9 1,5 4 35 4
8 140 100 23 26,2 0 26,2 1,5 8 35 4
9 142 102 23 27,35 0 27,4 1,5 8 35 4
10 142 92 23 28 0 28 1,5 8 35 4
11 145 102 23 29,95 0 30 1,5 8 35 4
12 198 177 23 32,2 0 32,2 1,5 8 35 4
Figure 3. Flow diagram of the developed software
Figure 4. A typical interface of developed software.
The developed software can also generating cross-sections throughout the 0Y axis where desired, and transform
them to orthogonal views using similar processes.
4 Case study and results
The highway rock cut at 110 km on the Konya-Antalya road in Turkey was chosen for the field experiment (Figure
5). In this region, the main formation is composed of limestone. Two primary discontinuity sets, which have
bedding and joints approximately perpendicular to each other, can be followed easily in the field with the naked
eye. The rock outcrop is relatively smooth, and discontinuities appear as fracture traces, not fracture surfaces.
The outcrop on which the data are recorded is shown in Figure 5(a), and the YZ plane of the experimental prism
and approximate discontinuity traces is shown in Figure 5(b).
The scan-line height was taken as 1.5 m. The outcrop orientation was 23°. The dimensions of the sampled prism
dimensions were selected as (x, y, z) 8×35×8 m. Then, the information system table (Table 2) was given to the
software as input, and all discontinuity planes were generated by the software.
The software can obtain and display several features from the selected sample prism, including all the
discontinuity planes (Figure 6(a)), all discontinuity planes with intersection lines (Figure 6(b)), the individual
intersection lines as dip direction vectors (Figure 6(c)), cross-sections throughout the 0Y axis where desired
(cross-section spacing were selected as 5 m; Figure 6(d)) and orthogonal views of the cross-sections (Figure
Figure 5. a) Outcrop chosen for the modeling (Konya-Antalya road at km 110 in Turkey), b) Approximate
discontinuity traces and sampled prism boundaries on the outcrop.
Figure 6. Model outputs: a) all discontinuity planes, b) all discontinuity planes with intersection lines, c) the
individual intersection lines as dip direction vectors, d) cross-sections throughout the 0Y axis where desired
(cross-section spacing were selected as 5 m), and e) orthogonal view of the cross-sections (2D view).
In this study, a rock structure was taken as a sample rectangular prism. Discontinuities were analyzed as linear
features within this prism. In the field, simple measurement instruments, including a tape measure and a
compass-clinometer, were adequate to collect the data from the exposures of rock masses. The recorded data
were used in the software directly and results were obtained easily and quickly. The primary functions of the
developed software are to: a) simulate the discontinuity network, b) obtain the intersection lines between
discontinuities, c) generate cross-sections where desired, directly from the simulated rock structure, and d)
provide a user-friendly interface. In spite of the sophisticated 3D model, our goal was to obtain useful, practical,
and rapid results.
Based on the results, we believe that the rock mass was adequately modeled for many practical engineering and
pre-evaluation applications. Despite assuming that the discontinuities have infinite sizes, the developed model
appears to be an improvement over other sampling methods such as photogrammetry or window sampling, and
thus the more realistic results can be obtained. By developing the present model, it is possible to analyze rock
masses for detailed engineering applications such as slope stability, fluid flow analysis, and blasting design.
This study has been supported by the Scientific Projects of Selçuk University (Turkey).
Andersson J. 1984. A stochastic model of a fractured rock conditioned by measured information. Water Resour Res 20(1):79–
Andersson J., Dverstorp B. 1987. Conditional simulations of fluid flow in three-dimensional networks of discrete fractures. Water
Resour Res 23(10):1,876–1,886.
Cundall P. A. 1971. Measurement and analysis of accelerations in rock slopes. Ph. D. thesis. London, Imperial College. 137.
Dershowitz W.S., Einstein H.H. 1987 Three-dimensional flow modelling in jointed rock masses. In: Proc. of 6th Cong. ISRM,
Montreal, Canada. Herget and Vongpaisal (eds.), Vol, 1, 87–92.
Elsworth D. 1986. A hybrid boundary-element-finite element analysis procedure for fluid flow simulation in fractured rock
masses. Int J Numer & Anal Meth Geomech 10.569–584.
Endo H.K., Long J.C.S., Wilson C.K., Witherspoon P.A.1984. A model for investigating mechanical transport in fractured media.
Water Resour Res 20(10):1,390–1,400.
Hocking G. 1977. Numerical Modeling in Rock Mechanics. Imperial Collage, London.
ISRM. 1978. Suggested Methods for Description of Discontinuities in Rock Mass: Int.Jour. Rock Mech Min Sci Geomech Abs, v.
Long J.C.S., Remer J.S., Wilson C.R., Witherspoon P.A. 1982. Porous media equivalents for networks of discontinuous
fractures. Water Resour Res 18(3):645–658.
Long, J.C.S., Gilmour, P., Witherspoon, P.A. 1985. A model for steady fluid flow in random three dimensional networks of discshaped
fractures. Water Resoul Res 21(8), 1105–1115.
Jing L. 1988. Formulation of discontinuous deformation analysis (DDA)—an implicit discrete element model for block systems.
Int J Eng Geol 1998;49:371–81.
Ming, C., K. Peter. 2006. Visualization of rock mass classification systems. Geotechnical and Geological Engineering, v. 24, p.
Pande G., Beer G., Williams J.R. 1990. Numerical Modeling in Rock Mechanics, John Wiley and Sons, 1990.
Shi G.H. 1996. Discontinuous deformation analysis technical note. First international forum on discontinuous deformation
analysis, June 12–14. Berkeley.
Smith L., Schwartz F.W. 1984. An analysis of the influence of fracture geometry on mass transport in fractured media. Water
Resour Res 20(9):1,241–1,252.
Williams J.R., Hocking G., Mustoe G.G.W., 1985. The Theoretical Basis of the Discrete Element Method, NUMETA 1985, Nume
Meth of Engin, Theory and Appl, A.A. Balkema, Rotterdam, January 1985.
Zhong, D.H., Li, M.C., Song, L.G. and G. Wang. 2006. Enhanced NURBS modeling and visualization for large 3D
geoengineering applications: An example from the jinping first-level hydropower engineering project, China. Computers &
Geoscience v.32, p. 1270-1282.
A New Approach to Rapid 3D Mapping of Rock Mass Structure