AB STRA CT The experiences of developing engineering geological models in karst areas for designing and construction purposes prove the necessity of considering at least three basic submodels: sedimentological, structural-tectonic and the weath- ering one. The research presented here deals with very important and frequently neglected segments in each of the submodels. Therefore, particular attention should be directed to: better understanding of carbonate sediment deposi- tion, determination of environment and diagenetic processes, study of the 3D anisotropy of discontinuity frequency, and differentiation of weathering zones. The given data and examples elaborate and justify such an approach, which enables a more realistic detailed engineering model, more reliable evaluations of the engineering geological/geotech- nical parameters and real site conditions. Keywords: engineering geological model, karst, limestone, anisotropy, heterogeneity, discontinuities, weather- ing, Croatia  Specific aspects of engineering- geological models in Croatian karst  Davor Pollak, Dražen Navratil and Tomislav Novosel Department of Hydrogeology and Engineering Geology, Croatian Geological Survey, Sachsova 2, HR-10000 Zagreb, Croatia; (davor.pollak@hgi-cgs.hr; drazen.navratil@hgi-cgs.hr; tomislav.novosel@hgi-cgs.hr) doi: 10.4154/gc.2013.17 1. INTRODUCTION The engineering geological model of the construction proj- ect area is a very important segment of geotechnical design. Creation of an engineering geological model for geotechni- cal designs is usually based on experience, known data, var- ious field observations and measurements, and assumes common engineering geological, laboratory testing, geo- physical surveying and borehole researching. This paper elaborates on the main geological issues observed at many investigation sites, while developing engineering geological models for the main motorways in the Croatian karst, at the Rijeka-Zagreb and Zagreb-Split-Dubrovnik sections having together more than 400 km in length (Fig. 1). The paper does not deal with detailed elaboration of every input data for future models, but with the generically important and sometimes neglected emphases that are found in building the engineering models for carbonate rocks. It is also the intention to document how the authors came to deal with the great heterogeneity of carbonate rock material and the heterogeneity/anisotropy of carbonate rock mass engi- neering properties. Geologia CroaticaGeologia Croatica Figure 1: Location map showing the recent distribution of Adriatic Carbon- ate Platform (AdCP) deposits (according : VLAHOVIĆ et al., 2005) and inves- tigated areas along new motorways. Geologia Croatica 66/3 219–232 13 Figs. 3 Tabs. Zagreb 2013 Geologia Croatica 66/3Geologia Croatica 220 2. GEOLOGICAL SETTING The Croatian karst area encompasses approximately 29,356 km2, or nearly 50% of the land surface of the Republic of Croatia. Carbonate rocks compose the central, coastal and submarine Croatian territory (Fig. 1). These sediments mainly originate from the Adriatic carbonate platform (AdCP). The large carbonate platform in this area started to develop during the Upper Triassic. During the Early Juras- sic period, it disintegrated into several smaller platforms, one of which was the AdCP (VLAHOVIĆ et al., 2005). Carbon- ate sedimentation in relatively calm and shallow marine con- ditions at the AdCP lasted during the Mesozoic until the Late Cretaceous. Afterwards, carbonate sedimentation was pro- longed under somewhat different conditions until the Eo- cene. The great thickness of carbonate sediments, (3500 -5000 m in Croatian territory), proves the constant subsidence of the carbonate platform during deposition (HERAK, 1991). Other factors, including the variability of the depositional environment and emersion periods, also testify to the sig- nificant dynamics of the whole carbonate platform. Tectonic activity between the Cretaceous and the Paleogene initiated regional emersion and termination of the platform sedimen- tation. Therefore during the Eocene, carbonate rocks were deposited at the margins of the disintegrated platform (VLAHOVIĆ et al., 2005). In addition, the flysch basin had been formed and the Dinarides uplifted. The uplifting of the Dinarides reached its maximum between the Oligocene and Miocene. The final result of formation of the Dinarides was an area with strongly faulted and crushed zones, reverse and overlapped contacts (HERAK, 1986). Since the beginning of the neotectonic episode, the Adriatic region continued to sink while the Dinarides uplifted (PRELOGOVIĆ et al., 2003). Later, more intensive uplifting caused erosion of the Dinarides and deposition of thicker sediments in depressions. The Quaternary deposits either cover Neogene sediments or are deposited directly on older bedrock. Karst is an area of specific hydrogeology and morphol- ogy which is characterized by sinkholes, caves and karst poljes, predominantly generated by the chemical weathering of carbonate rocks. The area of Dinaric karst has been inves- tigated since the middle of the 19th century and it is consid- ered as a „locus tipicus”. Therefore, some local terms for specific karst features are accepted worldwide. Mature and extreme karst with a great depth of karstifi- cation processes in Croatia is actually the consequence of lithology, tectonics, and action of water and climate condi- tions during geological history (HERAK et al., 1969; BA- HUN, 1974; PRELOGOVIĆ, 1975; HERAK & BAHUN, 1979; HERAK, 1991; FRITZ, 1991; VLAHOVIĆ et al., 1999; VLAHOVIĆ et al.; 2005). According to up to date knowledge, karstification of carbonate rocks in the AdCP occurred during many emersion periods in the Mesozoic. However, more intensive karstification on a regional scale occurred during the Upper Cretaceous – Tertiary, Middle and Late Eocene and Oligocene periods (HERAK et al., 1969). The Pleistocene age is very important for Dinaric karst de- velopment because the tectonic and hydrogeological proc- esses enabled the formation of karst poljes, (HERAK et al., 1969; PAVIČIĆ & RENIĆ, 1992) which are the largest karst phenomena. 3. THE BASIC POSTULATIONS FOR DEVELOPING AN ENGINEERING GEOLOGICAL MODEL The Croatian karstic area is a region of great variability of rock properties and geological situations affecting engineer ed design. This variability is three-dimensional and it affects all dimensions of the project from one geologic model unit to another. Therefore, before developing a model in a carbon- ate rock mass, an engineering geologist should first differ- entiate the geological elements that control the engineering properties and character of the construction ground: lithostra- ti graphic units, rock material types, structural-tectonic blocks, discontinuity sets, karst morphological features, weathering intensity and weathering „zones”. Therefore, field mapping activity needs to identify the expected critical data types how such data will be observed and collected so that rock mate- rial and rock mass properties can be noted and processed separately for each of the aforementioned elements. When this plan is completed, the engineering geological report will relate strictly to the design needs of the proposed engineered construction. Upon differentiation of the aforementioned geological elements, the authors propose the development of three mod- els: a sedimentological model, structural model and weath- ering model. The final engineering geological model will be the synthesis of these models. 3.1. Sedimentological model In carbonate rocks, most of the parameters set out above are investigated during sedimentological investigations mean- ing that the sedimentological model must be developed be- fore the engineering model (POLLAK & BRAUN; 1998 POLLAK, 2002). It is important to first define the expected lithostrati- graphical units of the project site or route alignment. The geological column should contain the basic descriptions of the major relevant sedimentological properties including: mineral composition, lithology, bedding thickness, interbeds, regularity of boundaries, structure, texture and porosity. In most cases, similar lithostratigraphic units can be merged according to their relevant engineering geological properties and conditions in order to set up engineering geological units. The usual method is to make a reconnaissance inspec- tion of the project area or route and to take representative laboratory samples of each apparent rock-unit material. These samples are than subjected to „index” testing in order to gain the first approximation of the gross engineering char- acter of each apparently different rock mass. It has been found that limestone texture can be ideally classified according to DUNHAM (1962), enabling a unified and simple description of rock material, according to depo- sitional texture observed directly in the field. Depositional textures, furthermore, then lead to recognition of the parent Pollak et al.: Specific aspects of engineering-geological models in Croatian karst Geologia Croatica 221 depositional environment and of later diagenetic processes. The importance of these factors in an engineering-geological modelling is paramount, because most of these effects have degraded (weakened) the limestone rock material. Considering the sedimentary nature of the limestone, each single unit is scrutinised for evidence of carbonate micro- facies. Units are often composed of several irregular vertical or lateral alterations of different microfacies, which cannot be considered as separate units. Therefore we use microscopic determination to reveal valuable details about the rock mate- rial: mineral and microfossil composition, texture, porosity, cementation, diagenesis, weathering and secondary defects including micro-fissures and cracks. The basic limestone and dolomite micro-textural properties are defined according to FOLK, (1959). The limestone microfacies determination is further done according to standard microfacies types (SMF) (FLÜGEL, 1982). The authors have found that carbonate rock materials in Croatian karst are frequently structurally modi- fied both by intensive diagenetic processes, as well as having frequent numeric defects due to intensive tectonics. When this form of „damage” is noted, the SMF should be modified and some new microfacies types added. 3.2. Structural model The structural model now is used as a means of defining the presence of structural blocks, determined on the basis of separate discontinuity sets (of joints, as planar discontinui- ties) orientation and their spacing. The discontinuities are found in mutually rhombic sets of two attitudes, each set caused by a single palaeotectonic deformational phase dur- ing geological history. These are the main planar breaks of weakness that define each structural block (also called struc- tural domains). In the investigated area of carbonate rocks, four tectonic phases that conditioned the evolution of tectonic blocks were present (MATIČEC, 2009). The most significant tectonic ac- tivities took place in the Upper Jurassic period, when the phase of compression tectonics started. In the late Cretaceous period the structures with a NNE-SSW strike were formed (Laramian orogene). In late Eocene and Oligocene the final uplift of the Dinarides mountain range began, along with formation of structures with strike NW-SE (Pyrenean oro- gene). Consequently, the Middle and Upper Miocene neo- tectonic phase followed, when regional stress assumed a gen- eral orientation N-S (MATIČEC, 2009). The relationships between the representative unit block (structural domains and major rock blocks) and the spatial po- sition of the linear infrastructure constructions are considered here. The construction of new linear infrastructure in bedrock terrain often produces new, fresh rock faces in road cuts which are of prime consideration. This new planar face sometimes creates the fourth face of a new tetrahedral block of rock that can be acted upon by gravity, and thus becomes unstable at the highway cut, or in a more limited sense, at the face or crown of a tunnel, during the construction process. The intensity of rock mass fracturing has fundamental importance regarding its great potential negative instability influence on rock mass characteristics before, during and af- ter excavation. Therefore, evaluation of the discontinuity fre- quency has usually been performed during the field investi- gations for any important building structure or mining project. In many areas, it is possible to collect such data ex- clusively by drilling. However, the usage of a single value for the linear frequency (discontinuity number determined per the unit line length in the defined direction through the rock mass), is generally not enough to ensure the correspond- ing intensity of the fracture quantification for the needs of the engineering design. For the rock mass which contains a definite number of sets of parallel, planar and persistent dis- continuities, HUDSON & PRIEST (1983) showed that the discontinuity frequency could be reduced to a single vector, as a rock mass property. For an example of the rock mass which contains N discontinuity sets, they also showed that their frequency along the scan line of any orientation is given with the expression: 1 cos N i s i il l q = = ∑ (1) where qi represents the angle between the scan line and the normal of the i-set of the discontinuity set, li is the frequency of the i-set of the discontinuity along the normal of that set where li is defined as the reciprocal value of the mean value of the discontinuity spacing 1 /x l= , and ls represents the total frequency along the determined scan direction. From this method of processing the orientation and fre- quency of discontinuity sets determined by the investigation field works, it is possible to determine the maximum and minimum values of the linear frequencies, their directions in three dimensions, i.e. to define the plane of the greatest anisotropy or define the afore mentioned values for the di- rections which are characteristic and important for a partic- ular linear object. Usage of the RQD index as one of the engineering meas- ures for classification of the fractured rock masses is wide- spread because of the easy measuring way on the outcrops and borehole cores. The analytical expression for the mutual relationship of RQD and random joint spacing, i.e. frequen- cies, is presented by PRIEST & HUDSON (1981) using the negative exponential distribution. HARRISON (1999) sets up the analytical determination of the threshold value which he prefers to the arbitrary thres- hold value which is used in practice (0.1m), and which is not sensitive to the whole range of the discontinuity spacing, i.e. frequencies which are met in practice. The use of the opti- mal RQD index threshold value greatly expands the range of RQD values encountered in rock mass and as such increa- ses the fidelity and hence the discriminatory nature of RQD. On the basis of the exponential model, according to PRIEST & HUDSON (1981) and the values of minimum lmin and maximum lmax discontinuity frequency, HARRI- SON (1999) suggested the function for determining the op- timal threshold value by means of the given formula: max max min min 2*t h l l l l   = ⋅  −   (2) Geologia Croatica 66/3Geologia Croatica 222 The simple classification of the rock masses for the pur- poses of engineering works such as tunnel construction, ex- cavations, foundations etc., (which can be done by means of the RQD values), has been known and used practically for a long time. What’s more, RQD presents an important pa- rameter for the categorization of the rock masses developed according to BIENIAWSKI (1973, 1989) and BARTON et al. (1974). The connection of RQD with the fundamental characteristics of the rock structure, i.e. their discontinuity frequency resulted in its widespread practical application. In literature, different statistical models for evaluation of the rock mass quality, PRIEST &HUDSON (1976), KUL- HAWY (1978), SEN (1984), have been suggested. These models were used for connecting RQD with the mean discontinuity number or with the average length of the intact core parts. One such model, which is based on the us- age of the negative exponential distribution for the model- ling of discontinuity spacing, was presented by PRIEST & HUDSON (1976). This model has been shown to be rather limited in the geological application (SEN, 1984), with one of the main limiting factors being the basic property of the exponential distribution that its mean value is equal to the standard deviation. Namely, there are numerous situations in nature where such equality is not valid. Because the dis- tribution of the intact parts lengths cannot be limited to the exponential distribution, the need to investigate other distri- butions is imposed. PRIEST & HUDSON (1976) and SEN (1984) showed that the analytical curves for RQD index opposite to l ob- tained on the base of the probability density function (p.d.f.) of the negative exponential distribution stand in the diagram almost always higher than the measured data, i.e. they over- estimate the value of the RQD index. SEN (1984) noticed that the standard deviation of the intact parts length plays an additional important role in the rock mass quality classifica- tion, i.e. an increase in the standard deviation improves the rock mass quality. He therefore suggested the log-normal distribution as the more appropriate model for connecting the RQD index with the average number of discontinuities, where the expected RQD is expressed by the formula: ( ) 1( ) 100 1 2Lnx Lnx Ln tE RQD F l s s    = − − +       (3); in which sLnx = standard deviation of the logarithmic intact lengths, l = mean number of discontinuities, t = threshold value and F(•) = standardized normal probability distribu- tion function. In this model E(RQD) and l are implicitly connected because the parameter l originates from the Pois- son distribution. Figure 2, initially published by SEN (1984), shows the relationship between RQD and l for two distributions (where t = 0.1 m): – lognormal distribution with different values of s Lnx (two parameter distribution); – exponential distribution (single parameter distribu- tion). Lognormal distribution shows that RQD cannot be judged solely on the basis of mean intact length but, addi- tionally, the standard deviation of the intact length must be considered. This property could not be identified with the negative exponential distribution. The mutual relationship of Rock Quality Designation (RQD) and the discontinuity frequency in carbonate rocks based on analytical models reported in the literature has also been considered and analyzed. 3.3. Weathering model The existing and generally used weathering engineering classi fication of karst ground (FOOKES & HAWKINS, 1988; WALTHAM & FOOKES, 2003) is mainly based on the number and size of karstic features and hydrogeological characteristics of the area. This classification is adequate and sufficient on a regional scale. An engineering geological model at a detailed scale, e.g. 1:1000 or 1:5000 requires ad- ditional thorough investigations of the weathering intensity of the rock mass regarding overburden depth. This aspect of weathering in karst is frequently neglected despite its obvi- ous indicators and earlier studies (DEERE & PATTON, 1971; NOVOSEL et al., 1980). Unlike other, indissoluble rock materials, carbonate rock material in this area demonstrates the minor influence of weathering to its mechanical and physical properties (POL- LAK, 2007), due to rapid and effortless weathering progress along numerous discontinuities. Therefore, most of the in- vestigations of the weathering processes in karst should be directed toward discontinuity properties from the surface zone into the interior. 4. RESULTS Despite the continuous deposition of carbonate rocks in shal- low and low energy environments, anomalous deposition conditions were encountered through the investigated area, which were acted upon by later diagenesis, intensive tectonic deformation, and weathering processes. It is our observation Figure 2: RQD variation with l for Lognormal p.d.f. Pollak et al.: Specific aspects of engineering-geological models in Croatian karst Geologia Croatica 223 that even subtle differences in the depositional environment of individual carbonate sediments, led to the creation of ma- jor heterogeneities and anisotropy that affect the final engi- neering of the highways. The authors here bring forth some important model find- ings, as gained from their route characterization work in car- bonate rock masses (Fig. 1). 4.1. Sedimentological model Detailed hand specimen evaluation, thin-section analysis, and the basic UCS testing usually explains the very high range of rock material mechanical properties and need for subdivision of the units. For example, it is very important to determine uniaxial compressive strength (UCS) of both crys- talline and stromatolitic segments from the Upper Triassic dolomite – "Hauptdolomit" (T3 2,3). The same dolomite in the Gorski kotar region, is present along more than 5 km of the Zagreb to Rijeka motorway. Stromatolitic dolomite repre- sents about 10–20% of the dolomite mass, but its lower uni- axial strength and laminated texture disrupts the homogene- ity and isotropy of the rest of the rock mass having a higher uniaxial strength (Tab. 1). The difference between the me- chanical properties of the crystalline and stromatolitic dolo- mite is significant, especially in weathered rock material. In the area between the Maslenica settlement and the "Sveti Rok" tunnel (13 km of motorway), microscopic anal- ysis confirmed the assumptions on close association of car- bonate rock material strength and its microfacieses (Figs. 3, 4). Despite the fact that just several samples were analyzed, some trends were noticeable. Firstly, a very strong and com- plete degree of calcite cementation made for a good and tight fit between the coarse grains of Neocomian limestone brec- cias (1d), which is the main reason for its greater strength in comparison to the other Lower Cretaceous (Neocomian) limestones (1a, b, c). The uniaxial compressive strength of recrystallized (2a) and dolomitized (2d) segments of Upper Cretaceous (Cenomanian) limestones has been significantly reduced in comparison to its parent materials (2b). Promina Table 1: Uniaxial compressive strength of different microfacies segments of "Hauptdolomite" in Gorski kotar region (POLLAK & BRAUN, 1998). DOLOMITES uniaxial strength / MPa texture (num. of samples) mean range CRYSTALLINE (46) 85 39-131 STROMATOLITE (8) 43 26-60 Figure 3: Subdivision of engineering geological units according to microfacies determination with appropriate microphoto. Zagreb-Split motorway, sec- tion Tunnel St. Rok-Maslenica. Geologia Croatica 66/3Geologia Croatica 224 carbonate rocks (4) consist of poorly sorted conglomerates (4a) and uniformly graded calcarenites (4b). While consid- ering these deposits, the significantly lower uniaxial com- pressive strength of conglomerates should be noted (Fig. 4). As opposed to the calcarenites, the matrix of the conglomer- ates is of variable quality and locally even porous which is, besides the grain size, considered to be the main reason for the lower strength of the conglomerate. At this point we were rewarded with the fact that deter- mination of subunits was possible on the basis of detailed field notes and some laboratory index property testing. In order to bridge the gap to the engineering geological recom- mendations, some of the sedimentological properties re- quired consideration from the engineering perspective, e.g. the bedding thickness. There are numerous examples to support the thesis of the influence of depositional and diagenetic conditions on the rock material and rock mass engineering properties. One of the most obvious examples is the influence of bedding thickness on block sizes of the carbonate rock mass. Bedding is the consequence of changes in physical, chem- ical and biological conditions of sedimentary facies and dia- genetic processes (COLLINSON & THOMPSON, 1989). Therefore a layer is considered to be a geological body sepa- rated from the surrounding rock mass by changes in its origi- nal grain-size distribution, mineralogical composition, litho- logical proprieties, orientation or pattern of particles or simply, by the termination of sedimentation. From the sedimentolog- ical perspective, the layer isn’t always bounded with discon- tinuities while considering the sedimentary rocks in engineer- ing sense it is usually very important to measure exclusively the spacing (frequency) of bedding discontinuities. The importance of such an approach can be confirmed in many segments, but here it is elaborated through determi- nation of the average field measured joint block size and shape in carbonate rocks. Based on the data from 32 mega lithostratigraphic units across the AdCP (3760 measurements of bedding thickness and 678 locations of block size calcu- lations), the analysis confirms the assumption that the „en- gineering bedding thickness” in the investigated area greatly affects the spacing of the secondary sets of discontinuities and naturally also the block sizes (Fig. 5). The diagram also shows the average block shapes, which are dominantly slightly tabular in well bedded, or equidimensional in the poorly bedded carbonate rock mass. 4.2. Structural model As previously stated, more than half of the Croatian territory is composed of carbonate deposits, i.e. a great number of linear infrastructural objects have been built in carbonate deposits in the last 15 years. Here, data collected in field in- ve stigations for the needs of geotechnical design are analys ed, by means of the models mentioned in the previous section, with the aim of improving the definition of the en gineering geological model. The collected data, which were obtained during the core drilling and engineering geological mapping for the design needs of important tunnels in carbonate de- posits, were analyzed (the segment of representative input data are presented in Table 2). Several lithographic units in the carbonate deposits were processed, where the blocks can generally be divided into equidimensional and tabular ones regarding the shape. In Table 3, the calculated minimum and maximum discontinu- ity frequencies, as well as their directions and planes of max- imal anisotropy are presented. Table 3 also contains the cal- culated values of the discontinuity frequencies based on the surface data in the direction of the performed drilling hole and strike of the designed tunnel axis. It is obvious from the data that the equidimensional blocks have lower maximum anisotropy regarding the discontinuity frequencies in differ- ent directions than the tabular blocks, which can be expected with regard to the block geometry itself. The values of the discontinuity frequencies and other pa- rameters from Tables 2 and 3 are presented in the following items, for two characteristic block types (equidimensional and tabular) which most often appear in carbonate rocks. Figure 6a presents the orientations of the separate dis- continuity sets and the axis of the tunnel "Mala Kapela" for structural unit A – equidimensional block (Fig. 6c). In Fig- ure 6b the 3D frequency values in all possible directions for the determined sets and their associated mean spacing are presented, the plane of maximum anisotropy is distinguished, Figure 4: The ranges and averages of uniaxial compressive strength (UCS) for particular microfacies (19 samples in total). Zagreb-Split motorway, sec- tion Tunnel Sveti rok-Maslenica. Figure 5: Correlation of the average bedding discontinuity spacing and block size in the Croatian carbonate rocks from the Middle Triassic to the Oligocene (continuous line). Dashed line represents ideal, calculated equi- dimensional blocks. Pollak et al.: Specific aspects of engineering-geological models in Croatian karst Geologia Croatica 225 Table 2: Example of discontinuity orientation, strike of tunnel axis and frequency data for two characteristic block types analyzed for the design needs of tunnels in carbonate rocks. Structure Structural unit Number of sets Discontinuity orientation (dip direction / dip angle) Normals of discontinuity sets (trend / plunge) Strike of tunnel axis Mean discontinuity spacing (m) Mean frequency (m–1) Tunnel “Mala Kapela” A 3 148/37 050/77 306/67 328/53 230/13 126/23 3 – 183 0.78 1.05 0.94 1.28 0.95 1.06 B 3 170/23 210/75 160/77 350/67 030/15 340/13 28 – 208 0.78 1.05 0.94 1.28 0.95 1.06 E 3 179/22 025/70 092/80 359/68 205/20 272/10 41 – 221 0.58 0.98 1.01 1.74 1.02 0.99 I 3 270/13 200/70 292/69 090/77 020/20 112/21 41 – 221 1.17 1.08 1.32 0.85 0.93 0.76 K 3 190/15 259/73 318/78 010/75 079/17 138/12 17 – 197 1.17 1.08 1.32 0.85 0.93 0.76 Tunnel “Umac” A 3 208/48 018/49 107/88 028/42 198/41 287/2 88 – 268 0.56 0.69 0.64 1.79 1.46 1.56 Tunnel “Brinje” D 3 59/13 241/77 182/79 239/77 061/13 002/11 18 – 198 0.21 0.89 0.54 4.78 1.13 1.84 D 4 59/13 241/77 182/79 129/84 239/77 061/13 002/11 309/6 18 – 198 0.21 0.89 0.54 0.54 4.78 1.13 1.84 1.84 E 3 103/13 248/80 181/89 283/77 068/10 001/01 20 – 200 0.17 0.78 0.86 5.97 1.27 1.17 E 4 103/13 248/80 181/89 300/83 283/77 068/10 001/01 120/07 20 – 200 0.17 0.78 0.86 0.86 5.97 1.27 1.17 1.17 Table 3: Global minimal and maximal frequency with their orientation and plane of maximal anisotropy, frequencies in direction of drilling, tunnelling direction and optimal threshold values for two characteristic block types of tunnels designed in carbonate rocks. Structure Structural unit / number of sets Min. frequency (m–1) Max. frequen- cy (m–1) Orientation of min. frequency (trend / plunge) Orientation of max. frequency (trend / plunge) Plane of extreme frequencies (max. anisotropy) Calculated frequency in direction of drilling (m–1) Mean frequencies from boreholes (m–1) Calculated frequency in tunnelling direction (m–1) Calculated optimal RQD threshold value (m) Tunnel “Mala Kapela” A / 3 1.00 2.17 220/10 290/20 281/20 1.65 4.6 1.79 1.32 B / 3 0.72 2.80 250/5 0/35 334/38 1.66 5.5 2.00 1.31 E / 3 0.92 2.55 115/10 265/55 202/73 2.13 8 2.02 1.25 I / 3 0.65 1.91 285/15 65/50 8/66 1.42 4.1 1.17 1.71 K / 3 0.64 1.89 175/15 95/40 103/40 1.25 3.25 1.02 1.74 Tunnel “Umac“ A / 3 1.52 3.03 195/45 55/5 142/59 2.21 3.45 2.52 0.91 Tunnel “Brinje” D / 3 1.02 5.56 95/10 5/70 9/70 5.26 10.1 3.35 0.75 D / 4 2.48 6.33 95/10 335/60 10/65 5.46 10.1 4.01 0.48 E / 3 1.27 6.36 145/10 355/70 57/80 6.06 6.4 2.11 0.63 E / 4 2.15 6.66 30/5 295/55 304/55 6.20 6.4 2.27 0.50 Geologia Croatica 66/3Geologia Croatica 226 and the direction of borehole as well as the strike of the tun- nel axis are marked. Figure 7 provides a similar presentation for the tabular block (Fig. 7c), i.e. structural unit E in the "Brinje" tunnel (Fig. 7a and b). Figure 7b), which evidently shows a greater anisotropy of frequency for the tabular blocks. According to the equation [2], for the calculation of an optimal threshold value of RQD, it is necessary to determine the value of the minimum and the maximum discontinuity frequency (l min and l max), with regard to the given orien- tations of the defined sets and their spacing (Table 3). Fig- ures 8 and 9 present the polar plot of discontinuity frequency in the plane of extreme frequencies and differences in the discriminatory nature of the RQD index for the customary thres hold value and calculated optimal threshold values for the structural units analysed in figures 6 and 7. Figure 6: The "Mala Kapela" Tunnel – structural unit A (3 sets): a) Equal angle lower hemispherical projection of example data, b) three-dimensional locus of discontinuity frequency with plane of extreme (range) frequencies, c) equidimensional blocks. Figure 7: The "Brinje" Tunnel – structural unit E (3 sets): a) Equal angle lower hemispherical projection of example data, b) three-dimensional locus of discontinuity frequency with the plane of extreme frequencies, c) tabular blocks. Pollak et al.: Specific aspects of engineering-geological models in Croatian karst Geologia Croatica 227 To establish the connection between the RDQ index and the average number of discontinuities (frequency), the afore mentioned exponential and log-normal model was applied. The data was collected from cores of 46 boreholes, drilled during the design phase of an important tunnel, through mas- sive to thick layered limestones of Cretaceous age. The rea- son those rocks were chosen is because the data on all the intact length (not only of longer than 0.1 m) were available, which is necessary for the estimation of s Lnx. Highly frac- tured and cavernous zones were excluded during the calcu- lation of RQD and average number of discontinuities. Results of the analysis confirm the justification of the log-normal model usage, as opposed to the exponential one, for the purpose of connecting these two parameters (RQD vs. average number of discontinuities – l) (Fig. 10). Log- normal function of the probability density is defined with the formula [3] where the standard deviation of the logarithmic intact lengths s Lnx = 0.78. 4.3. Weathering model Despite the enormous importance of karstification processes on engineering problems in the Croatian karst, that segment of engineering-geological model is not the focus of this paper. The main superficial (sinkholes) and subterranean (ca ves, caverns) karstic features have been studied by many authors (ŠARIN & HRELIĆ, 1984; KAPELJ, 1985; BOŽIČEVIĆ et al., 1995; JAMIČIĆ & NOVOSEL, 1999; WALTHAM & FOOKES, 2003; KAPELJ et al., 2004; GARAŠIĆ, 2005; Figure 8: The "Mala Kapela" Tunnel – structural unit A (3 sets): a) Polar plot of discontinuity frequency; b) polar plot of RQD using customary threshold value 0.1 m and optimal threshold value. Figure 9: The "Brinje" Tunnel – structural unit E (3 sets): a) Polar plot of discontinuity frequency, b) polar plot of RQD using customary threshold value 0.1 m and optimal threshold value. Geologia Croatica 66/3Geologia Croatica 228 BOSTJANČIĆ et al., 2011) and their results should be ac- knowledged in engineering design encompassing karstified media. Here, the influence of weathering processes on gen- eral discontinuitiy properties and rock mass quality is con- sidered. The carbonate rocks in Croatia are predominantly very "pure", with more than 92% carbonate mineral content. Car- bonate composition and other properties of the area, includ- ing low primary porosity of rock material (< 5%), the thick- ness of the carbonate sediment complex, tectonics and climatic conditions enabled intensive weathering of the car- bonate rocks and the formation of karst. The probable "zoning" of the carbonate rock mass should be examined during the investigation phase for the engineering geological model. Initially that is usually done by thorough engineering geological mapping. In the later phases the "zoning" should be verified with geophysical in- vestigations and core drilling. Most of the analyzed geophysical profiles which are done in the investigation phase for the Croatian motorways across the AdCP display surprisingly ‘regular’ zones. The compilation histogram of 54 shallow refraction seismic pro- files with total length of 4160 m, displays clear zones in the limestone rock mass (Fig. 11). Regarding the relatively shal- low maximum reach of several decametres in most analyzed geophysical profiles, three zones can be noted: SWZ – sur- face weathering zone; UWZ – upper weathering zone; LWZ – lower weathering zone. It is also obvious that fresh rock mass, with expected Vp > 6000 m/s is rarely reached as a minor segment in LWZ. Despite the clearly visible zoning in Figure 11, determi- nation of the weathering zones at the specific location shouldn’t be based on geophysical investigations alone. The weathering processes and the discontinuity properties char- acteristic for each weathering zone can usually be recognized in the borehole core. Besides the obvious variation of weathering marks at the discontinuities in borehole core, it is also noted that RQD values can show a trend of block size change as an indica- tion of weathering intensity related to the rock mass depth. From the extensive amount of borehole data (based on 2930 intervals in 463 boreholes with the total drilled interval of 4877 m), the trend of block size increase with depth is clearly visible (Fig. 12). Still, it is important to say that such clear zoning is rarely visible in boreholes on smaller and extremely karstified areas. The ‘zoning’ of engineering properties in the karst area is mainly a consequence of changes in the width and infill- ing of discontinuities. Namely, from the fresh rock to the strongly weathered rock mass, the width of discontinuities increases significantly but the infilling also generally changes from none, over carbonate through clayey to none. This leads to the derivation of characteristic weathering zones in the Croatian karst, which can usually be incorporated in the en- gineering-geological model (Fig. 13). The definition of the weathering zones in this model are mainly based on these two factors – the width and infilling of discontinuities. Figure 10: Variation of the RQD index versus l for Exponential and Lognor- mal probability density function (p.d.f.) in Croatian carbonate rocks (Upper Juras sic up to Upper Cretaceous limestones and dolomites). Figure 11: The distribution of the primary seismic waves velocities (Vp) in different weathering zones in limestones across the AdCP. Figure 12: The histogram of RQD index values in relation to the depth of the drilled interval from boreholes across AdCP. Pollak et al.: Specific aspects of engineering-geological models in Croatian karst Geologia Croatica 229 5. DISCUSSION The main objective during the engineering geological pro- cessing of the field and laboratory data is the synthesis of geological data and quantitative parameters for use by the design engineers. The synthesis of all the available data should produce an engineering geological model, which con- tains all the important engineering data and presents it in simple and practical way. The authors advocate that beside the studies of the char- acteristic karstic morphological features, there are other very important segments of the engineering geological model of karstified rocks: the sedimentological, structural and weather- ing model; which enables a better and more reliable definition of weathered rock mass properties and quality. In that sense, the engineering geologist’s eye must seek and appreciate min- eral composition, formation, diagenesis, rock material texture, structural discontinuities and the evidence of alteration and/ or weathering which tends to degrade the basic engineering properties and parameters of each rock-mass unit. Engineering sedimentological modelling in carbonate rocks begins with detailed field observations and micro- scopic investigations of samples selected as being represent- ative of each apparent three-dimensional rock mass unit. Here the field-data collection needs are: unit lithology, geometry of the sediment body, texture, structure, bedding plane properties and bedding thickness. At the beginning of the field mapping and data collection stage, there is the nat- ural challenge of making a preliminary evaluation of the likely impact of different factors on the engineering charac- ter of rock masses to be encountered. Most of the carbonate rock mass in Croatia is bedded and the rock mass therefore has a basic anisotropy, and this requires constant attention to the need to identify three-di- mensional block parameters, particularly as it affects the di- rection not only of motorways, but of all other linear infra- structural projects. With regard to the anisotropy of the discontinuity frequency in different directions, the question is whether the vertical boreholes drilled during the geotech- Figure 13: Characteristic weathering zones in Croa tian karst with corresponding basic descrip- tions. (Photo examples are from the "Mala Kapela" tunnel area). Geologia Croatica 66/3Geologia Croatica 230 nical design of the tunnels give a real evaluation of the rock mass, i.e. if they underestimate or overestimate the fractur- ing of the rock mass in the direction of tunnel excavation? According to the presented data, the anisotropy of discon- tinuity frequency in equidimensional blocks is insignificant while in tabular blocks it is strongly expressed, especially in the case of the usual three sets of discontinuities, each result- ing from a separate but severe historical phase of regional tec- tonism. Taking into consideration the realistic example (Fig. 6), the vertical core drilled at the "Brinje" tunnel underesti- mates the rock mass quality in the direction of the tunnel ex- cavation concerning the fracturing state. In addition to the calculation, this was proven by comparing core logging data with data obtained during excavation of the tunnel. Also, the mean values of the discontinuity frequencies obtained by core logging, in the majority of examples pre- sented in Table 3; do not match the calculated frequencies in the drilling direction. The reasons for this can be manifold: 1) the "in situ" rock mass is disturbed during the drilling process, i.e. the obtained core contains more fractures than are visible on the walls of the borehole due to fracturing or opening along cracks where weakened part in the intact rock exists; 2) the objective impossibility of collecting enough qualitative data concerning the discontinuity orientations and their spacing on the surface. Therefore, it is not possible to reliably define the mean values of the set discontinuity ori- entations, i.e. the size and shape of the unit block, only on the basis of the surface data. The impossibility of collecting sufficient qualitative data is often connected with weather- ing in the surface zones, coverage by Quaternary deposits and with a smaller blocks size on the surface, caused by the process of chemical weathering of the carbonate rocks. Er- rors in defining the aforementioned parameters directly in- fluence the results of the calculated frequencies obtained in the determined direction. So the most appropriate model should consist of reliably defined mean values of sets orientations and their frequencies. Therefore, it is necessary to bear in mind that the model for the analysis of the discontinuity frequency is based on three basic assumptions: • rock mass structure is composed of planar sets of par- allel discontinuities; • all discontinuities are strictly planar, which is not the true case; • each discontinuity is persistent along the observed area. Using the analysis of the optimum threshold value of the RQD index, according to HARRISON (1999) it is clear that its usage maximises the range of values of the RQD index, which is marked in the rock mass and as such it increases the reliability and consequently the discriminative nature of the RQD index. It is also clear from the diagram (Fig. 10) that the expo- nential model overestimates the RQD values for the analysed data, while the log-normal model is more justifiable regard- ing the security in geotechnical design in carbonate rock masses. However, it is necessary to have all the intact lengths for the estimation of s Lnx. Given the specific weathering mechanisms, the weath- ering zones in extremely karstified areas are frequently ir- regular. Therefore, the expression "weathering zone" shouldn’t be literally used in karstic regions. It is often the case that the separated environments mutually interlace and irregu- larly exchange vertically and laterally. Irregular "zoning" appears on the large faults or in tectonically fractured areas, where borders with other zones could be vertical or even in- verse. Depending on the geological properties of the inves- tigated area, the weathering zones in karst have very differ- ent characteristics, spreading and mutual relationships. Therefore, it is difficult to judge the differentiation of the "zones" of variable degrees of weathering in the carbon- ate rock mass, because of its great irregularity in every sense. There are also usually great differences in discontinuity prop- erties from the outcrop zone to the greater overburden depth of several hundred metres which should be quantified in the engineering geological model. In general, weathering is most intensive at the surface zone where mainly the chemical but also physical action am- plifies the appearance of the joint sets, which then appear to be open, and there are usually frequent wide discontinuities. In the areas with deeper overburden, the same discontinuity sets can be closed or even disappear. Therefore, in the car- bonate rock mass one should expect an increase of block size (and a related decrease in spacing frequency) with the rise of overburden depth. This is visible in numerous cuttings along the newly built motorways in Croatia. 6. CONCLUSIONS Beside the general and well known procedures for building the engineering geological model in any rock mass, a car- bonate rock mass demands some additional perspectives. If we adopt the regional model which quantifies the karstifica- tion intensity of the area, the detailed investigations should also enable development of three main submodels: 1. sedi- mentological, 2. structural and 3. weathering. Here, some very important aspects of each submodel are studied and presented respectively. 1. The total understanding of the history of sedimenta- tion provides explanation and answers about engineering behaviour and properties of both rock material and rock mass. Therefore, the sedimentological investigations of car- bonate rocks are greatly improving the creation of the engi- neering geological model because they enable more precise and logical differentiation of engineering units through the described procedure. The microfacies are characteristic for strictly defined depositional environments and later geolog- ical conditions and they are proven to have a major influence on the many engineering properties of carbonate rock mate- rial and mass. Therefore, the definition of microfacies usu- ally leads to better comprehension of the variability in me- chanical properties in carbonate rock material of each of the engineering geological units. 2. The determination of the extreme values of the dis- continuity frequencies has successfully been applied in plan- ning the sampling strategy by drilling, different classifica- Pollak et al.: Specific aspects of engineering-geological models in Croatian karst Geologia Croatica 231 tions and general understanding of the mechanical and hydrological characteristics of the fractured carbonate rock masses. Based on the considerable quantity of data from the boreholes in the massive to thick layered limestones of Cre- taceous age, the justification supposition of the log-normal model usage was confirmed as opposed to the exponential one with the purpose of RQD connecting with the average number of the discontinuities. 3. Despite clearly visible weathering zoning in the car- bonate rock mass on a great amount of data it is very impor- tant to stress that weathering "zones" in Croatian karst are very irregular in space and in almost every other parameter, as it is typical for all Dinaric karst areas. Still, creation of the detailed and valid engineering geological model of a karsti- fied area should include consideration of possible weather- ing zones and its engineering properties. Four weathering zones regarding the properties of discontinuities and weath- ering marks are characteristic: fresh rock mass (I, FRM), lower weathering zone (II, LWZ), upper weathering zone (III, UWZ) and surface weathering zone (IV, SWZ). Determination of the geological genesis of each three- dimensional rock mass is the key to understanding the over- all condition of the rock for present construction. Once com- pleted, this leads to the most accurate determination of the rock material and rock mass engineering properties and pre- diction of post-engineering behaviour. REFERENCES BAHUN, S. (1974): Tektogeneza Velebita i postanak Jelar-naslaga [Tec- togenesis of the Velebit and genesis of Jelar deposits – in Croatian]. – Geol. vjesnik, 27, 35–51. BARTON, N., LIEN, R. & LUNDE, J. (1974): Engineering classifica- tion of rock masses for the design of tunnel support.– Rock Me- chanics, 6, 183–236. BIENIAWSKI, Z.T. (1973): Engineering classification of jointed rock masses.– Transactions of the South African Institution of Civil En- gineers, 15, 335–344. BIENIAWSKI, Z.T. (1989): Engineering rock mass classifications.– John Wiley & Sons, New York, 251 p. BOSTJANČIĆ, I., POLLAK, D., PODOLSZKI, L. & GULAM, V. (2011): The geological setting of the sinkholes in bare Croatian karst.– In: LAVEROV, N.P. & OSIPOV, V.I. (eds.): EngeoPro, En- vironmental Geosciences and Engineering Survey for Territory Pro- tection and Population Safety, Proceedings, Moscow, 179–184. BOŽIČEVIĆ, S., ŠABAN, B. & VULIĆ, Ž. (1995): Podzemne šupljine u temeljima objekata na području Rijeke [Underground cavities in the foundations of buildings in Rijeka – in Croatian]. – In: VLAHOVIĆ, I., VELIĆ, I., ŠPARICA, M. (eds.): First Croatian Geological Congress, Proceedings, 1, 121–124. COLLINSON, J.D. & THOMPSON, D.B. (1989): Sedimentary struc- tures.– Unwin Hyman, Boston-Sydney-Wellington, 207 p. DEERE, D.U. & PATTON, F.D. (1971): Slope stability in residual soils.– 4th Pan-American Conference on Soil Mechanics and Foundation Engineering, Proceedings, San Juan, Puerto Rico, 87–170. DUNHAM, R.J. (1962): Classification of the carbonate rocks according to depositional texture – Classification of carbonate rocks American association of petroleum geologists Memoirs, 1, 108–121. FLÜGEL, E. (1982): Microfacies analysis of limestones.– Springer-Ver- lag, Berlin-Heidelberg, 633 p. FOLK, R.L. (1959): Practical petrographic classification of limestones – Bulletin of American Association of Petroleum Geologysts, 43, 1–38. FOOKES, P.G. & HAWKINS, A.B. (1988): Limestone weathering: its engineering significance and a proposed classification scheme.– Quarterly Journal of Engineering Geology, 21, 7–31. FRITZ, F. (1991): Utjecaj recentnog okršavanja na zahvaćanje voda [The impact of recent karstification on water extraction – in Croatian]. – Geol. vjesnik, 44, 281–288. GARAŠIĆ, M. (2005): Speleological research of caverns in trace of the highways in Croatian karst.– In: VELIĆ, I., VLAHOVIĆ, I., BIONDIĆ, R. (eds.): Third Croatian Geological Congress, Ab- stracts book, Opatija, 191–192. HARRISON, J.P. (1999): Selection of threshold value in RQD assess- ments.– International Journal of Rock Mechanics and Mining Sci- ences, 36, 673–685. HERAK, M. (1986): A new concept of geotectonics of the Dinarides.– Acta Geologica, Croatian Academy of Sciences and Arts, 21, 35– 83. HERAK, M. (1991): Dinaridi: Mobilistički osvrt na genezu i strukturu [Dinarides: Mobilistic review of the genesis and structure – in Croatian].– Acta geologica, 21/2, 35–117. HERAK, M. & BAHUN, S. (1979): The role of the calcareous breccias (jelar Formation) in the tectonic interpretation of the High Karst Zone of the Dinarides.– Geol. vjesnik, 31, 49–59. HERAK, M., BAHUN, S. & MAGDALENIĆ, A. (1969): Pozitivni i negativni utjecaji na razvoj krša u Hrvatskoj [Positive and negative influences on the development of karst in Croatia – in Croatian].– Krš Jugoslavije, 6, 45–78. HUDSON, J.A. & PRIEST, S.D. (1983): Discontinuity Frequency in Rock Masses.– International Journal of Rock Mechanics and Mi- neral Science & Geomechanics Abstracts, 20/2, 73–89. JAMIČIĆ, D. & NOVOSEL, T. (1999): The dynamics of tectonic mod- eling of some caves in the karst region (Croatia).– Geol. Croat., 52/2, 197–202. KAPELJ, J. (1985): Neki rezultati primjene geomorfološke analize relje- fa u uvjetima krškog terena [Some results of the application of geo- morphic analysis of relief in karst terrain – in Croatian].– Geol. vjesnik, 38, 121–129. KAPELJ, J., KAPELJ, S. & SINGER, D. (2004): Spatial distribution of dolinas and its significance for groundwater protection in karst ter- rain.– XXXIII International Association of Hydrogeologists Con- gress, Proceedings, Mexico, 1–4. KULHAWY, F.H. (1978): Geomechanical Model for Rock Foundation Settlement.– Journal of Geotehnical Engineering Division, ASCE, 104, 211-228. MATIČEC, D. (2009): Krški Dinaridi [Karst Dinarides – in Croatian] – In: VELIĆ, I. & VLAHOVIĆ, I. (eds.): Tumač Geološke karte Republike Hrvatske 1:300000, Hrvatski geološki institut, Zagreb, 105–106 p. NOVOSEL, T., TUŠAR, Z., MULABDIĆ, M., GARAŠIĆ, M. & KORAŽIJA, S. (1980): Ocjena stabilnosti kosina u zasjecima (us- jecima) građenih od karbonatnih stijena [Rating stability in notches (cuttings) built of carbonate rocks – in Croatian].– Zbornik radova 5. Simpozija jugoslavenskog društva za mehaniku stijena i podzemne radove, 1, 185–192. PAVIČIĆ, A. & RENIĆ, A. (1992): Recent karstification and water drain- age from Bakovac Valley (Lika, Croatia).– Geomorphology and Sea, Proceedings of the International Symposium, Mali Lošinj, 127–133. POLLAK, D. (2002): Ovisnost inženjerskogeoloških svojstava karbon- atnih stijena o njihovim sedimento-petrološkim značajkama (Trasa Jadranske autoceste: "Tunel Sveti Rok – Maslenica") [Dependence of engineering-geological properties of carbonate rocks on their Geologia Croatica 66/3Geologia Croatica 232 sedimentary-petrological characteristics (Adriatic highway sec- tion: "Tunnel Sv. Rok – Maslenica") – in Croatian].– Unpubl. MSc Thesis, Faculty of Mining, Geology and Petroleum Engineering, Zagreb University, 108 p. POLLAK, D. (2007): Utjecaj trošenja karbonatnih stijenskih masa na njihova inženjerskogeološka svojstva [Influence of the Carbonate Rock Masses Weathering on Their Engineering-Geological Prop- erties – in Croatian].– Unpubl. PhD Thesis, Faculty of Mining, Ge- ology and Petroleum Engineering, Zagreb University, 299 p. POLLAK, D. & BRAUN, K. (1998): Sedimentology in the service of en- gineering geology: Study of some results of the explorations for high- way construction and tunneling in Croatia.– In: MOORE, D. & HUNGR, O. (ed.): Proceedings, Rotterdam: A. A. Balkema, 195–199. PRELOGOVIĆ, E. (1975): Neotektonska karta SR Hrvatske [Neotec- tonic map of the Croatia – in Croatian].– Geol. vjesnik, 28, 97–108. PRELOGOVIĆ, E., PRIBIČEVIĆ, B., IVKOVIĆ, Ž., DRAGIČEVIĆ, I., BULJAN, R. & TOMLJENOVIĆ, B. (2003): Recent structural fabric of the Dinarides and tectonically active zones important for petroleum-geological exploration in Croatia.– Nafta, 55/4, 155– 161. PRIEST, S.D. & HUDSON, J.A. (1976): Discontinuity Spacings in Rock – International Journal of Rock Mechanics and Mineral Science & Geomechanics Abstracts, 13, 135–148. PRIEST, S.D. & HUDSON, J.A. (1981): Estimation of Discontinuity Spacings and Trace Length Using Scan-line Survey.– International Journal of Rock Mechanics and Mining Sciences, 18, 183–197. ŠARIN, A. & HRELIĆ, Đ. (1984): Ocjena hidrogeoloških svojstava kar- bonatnih naslaga na temelju gustoće ponikava u niskom i dijelu vi- sokog krša Hrvatske [The evaluation of hydrogeological properties of carbonate sediments on the basis of sinkhole density in low and high karst of Croatia – in Croatian].– VII Jugoslavenski simpozij, Budva, Zbornik radova, Knjiga 1, 117–125. SEN, Z. (1984): RQD models and fracture spacing – Journal of Geotech- nical Engineering, ASCE 110/2, 203–216. VLAHOVIĆ, I., VELIĆ, I., TIŠLJAR, J. & MATIČEC, D. (1999): Li- thology and Origin of Tertiary Jelar Breccia within the Framework of Tectogenesis of Dinarides - Some carbonate and clastic Succes- sions of the External Dinarides.– Harlod Reading’s IAS Lecture Tour, 23–25. VLAHOVIĆ, I., TIŠLJAR, J., VELIĆ, I. & MATIČEC, D. (2005): Evo- lution of the Adriatic Carbonate Platform: Palaeography, main events and depositional dynamics.– Palaeography, Palaeoclimatol- ogy, Palaeoecology, 220, 333–360. WALTHAM, A.C. & FOOKES, P.G. (2003): Engineering classification of karst ground conditions.– Quarterly Journal of Engineering Ge- ology and Hydrogeology, 36, 101–118. Manuscript received February 08, 2013 Revised manuscript accepted September 26, 2013 Available online October 28, 2013