с Г ISSN 1854-0171 / z £ «г # 4# У # # fjpjpj* J м/ л # # v-О ACTA GEOTECHNICA SLOVENICA ISSN: 1854-0171 Ustanovitelji Founders Univerza v Mariboru, Fakulteta za gradbeništvo, prometno inženirstvo in arhitekturo University of Maribor, Faculty of Civil Engineering, Transportation Engineering and Architecture Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo University of Ljubljana, Faculty of Civil and Geodetic Engineering Univerza v Ljubljani, Naravoslovnotehniška fakulteta University of Ljubljana, Faculty of Natural Sciences and Engineering Slovensko geotehniško društvo Slovenian Geotechnical Society Društvo za podzemne in geotehniške konstrukcije Society for Underground and Geotechnical Constructions Izdajatelj Publisher Univerza v Mariboru, Fakulteta za gradbeništvo, prometno inženirstvo in arhitekturo Faculty of Civil Engineering, Transportation Engineering and Architecture Odgovorni urednik Editor-in-Chief Bojana Dolinar University of Maribor Uredniki Co-Editors Jakob Likar Janko Logar Borut Macuh Stanislav Škrabl Bojan Žlender University of Ljubljana University of Ljubljana University of Maribor University of Maribor University of Maribor Posvetovalni uredniki Advisory Editors Heinz Brandl Chandrakant. S. Desai Pedro Seco e Pinto Vienna University of Technology University of Arizona National Laboratory of Civil Eng. Lektor Proof-Reader Paul McGuiness Naklada Circulation 250 izvodov - issues Cena Price 25 EUR/letnik - 25 EUR/vol.; (50 EUR for institutions/za institucije) Tisk Print Tiskarna Saje Revija redno izhaja dvakrat letno. Članki v reviji so recenzirani s strani priznanih mednarodnih strokovnjakov. Baze podatkov v katerih je revija indeksirana: SCIE - Science Citation Index Expanded, JCR - Journal Citation Reports / Science Edition, ICONDA - The international Construction database, GeoRef. Izid publikacije je finančno podprla Javna agencija za raziskovalno dejavnost Republike Slovenije iz naslova razpisa za sofinanciranje domačih periodičnih publikacij. Uredniški odbor Amin Barari Theodoros Hatzigogos Vojkan Jovičič Rolf Katzenbach Nasser Khalili Bojan Majes Svetlana Melentijevic Ana Petkovšek Borut Petkovšek Mihael Ribičič César Sagaseta Patrick Selvadurai Stephan Semprich Devendra Narain Singh Abdul-Hamid Soubra Kiichi Suzuki Antun Szavits-Nossan Kosta Urumović Ivan Vaniček Editorial Board Aalborg University Aristotle University of Thessaloniki IRGO-Ljubljana, President of the SloGeD Technical University Darmstadt The University of New South Wales, Sydney University of Ljubljana RodioKronsa, Soletanche-Bachy Group, Spain University of Ljubljana Slovenian National Building and Civil Engineering Institute University of Ljubljana University of Cantabria McGill University University of Technology Graz Indian Institute of Technology, Bombay University of Nantes Saitama University University of Zagreb Croatian geological survey Czech Technical University in Prague Naslov uredništva Address ACTA GEOTECHNICA SLOVENICA Univerza v Mariboru, Fakulteta za gradbeništvo, prometno inženirstvo in arhitekturo Smetanova ulica 17, 2000 Maribor, Slovenija Telefon / Telephone: +386 (0)2 22 94 300 Faks / Fax: +386 (0)2 25 24 179 E-pošta / E-mail: ags@uni-mb.si Spletni naslov web Address http://www. fg.uni-mb.si/journ al-ags/ The journal is published twice a year. Papers are peer reviewed by renowned international experts. Indexation data bases of the journal: SCIE - Science Citation Index Expanded, JCR - Journal Citation Reports / Science Edition, ICONDA- The international Construction database, GeoRef. The publication was financially supported by Slovenian Research Agency according to the Tender for co-financing of domestic periodicals. Uvodnik Editorial T. Fabjan in drugi Numerično modeliranje mehanskega obnašanja intaktne kamnine z modelom kontinuuma in z modelom z Voronoijevi-mi elementi - občutljivostna analiza T. Fabjan et al. Numerical simulation of intact rock behaviour via the continuum and Vor-onoi tessellation models: a sensitivity к analysis 4 ▲ 24 I. Ružić in drugi Ocena stabilnosti obalnih pečin z uporabo digitalnih posnetkov I. Ružić et al. A stability assessment of coastal cliffs using digital imagery F. Jiang in drugi Nova metoda za preizkušanje trdnosti gline na porušitev zaradi precejanja pri visokih vodnih tlakih F. Jiang et al. A new method for testing the anti-permeability strength of clay failure under a high water pressure 36 A 46 M. M. Hajitaheri & M. Hassanlourad Numerično modeliranje negativnega trenja plašča na posamičnem vertikalnem in nagnjenem pilotu M. M. Hajitaheri & M. Hassanlourad Numerical modeling of the negative skin friction on single vertical and batter pile Y. Zhuang & X. Cui Y. Zhuang & X. Cui Nasip na mehkih tleh ojačan s piloti za železnice visokih hitrosti: numerična in analitična raziskava Reinforced piled embankment for a high-speed railway over soft soil: a ^umerica! and analytical investigation 56 A J. C. Leal in drugi J. C. Leal et al. ■■■ Suha frakcija nezasičenih tal The dry fraction of unsaturated soils 66 P. Han in drugi P. Han et al. 76 Korozijski mehanizmi cementiranih zemljin za tri različne raztopine sulfatov Corrosion mechanisms for cemented soils in three different sulfate solutions ▲ Navodila avtorjem Instructions for authors 86 UVODNIK Z veseljem ugotavljamo, da se je v tem letu število prispelih člankov v uredništvo revije Acta Geotechnica Slovenica zelo povečalo, kar kaže na večjo razpoznavnost revije in tudi cenjenost med avtorji prispevkov. V uredništvu smo veseli takega odziva avtorjev, saj nam to omogoča izbor najboljših člankov. Istočasno pa moramo zavrniti vrsto zelo kvalitetnih prispevkov, ki bi bili vredni objave, vendar nam sedanji obseg revije tega ne dopušča. Izberemo lahko le približno deset odstotkov prispelih del, ostalim avtorjem dobrih člankov pa predlagamo objavo v drugih revijah. Medtem ko smo v lanskem letu že povečali število objavljenih člankov z deset na dvanajst, smo v letu 2015 objavili trinajst člankov. S tem povečujemo tudi stroške izdajanja revije, kar pa za Univerzo v Mariboru, Fakulteto za gradbeništvo, prometno inženirstvo in arhitekturo, kot izdajatelje, predstavlja velik problem, posebno še, ker se je sofinanciranje revije s strani Agencije za raziskovalno dejavnost Republike Slovenije v tem letu drastično zmanjšalo. Kljub temu pa ne želimo del stroškov izdajanja revije prenesti na avtorje, saj bi jih s tem postavili v neenak položaj glede možnosti objave člankov. V lanskem letu smo se v uredništvu odločili za sprejem novih članov v Uredniški odbor revije. Pri izboru upoštevamo predvsem kvaliteto dosedanjega dela posameznikov pri recenziranju člankov in njihovo pripravljenost za aktivno sodelovanje tudi v bodoče. Tako smo v lanskem letu vključili dr. Ano Petkovšek, ki prihaja s Fakultete za gradbeništvo in geodezijo v Ljubljani, Katedre za mehaniko tal z laboratorijem in ima izkušnje predvsem s področja geomehanike, gradnje prometnic in geotehnike okolja. V tem letu se nam je pridružil dr. Kosta Urumović s Hrvaškega geološkega inštituta, Zavoda za hidrogeologijo in inženirsko geologijo. Njegovo delo na HGI je osredotočeno predvsem na področje hidrogeologije, na katerem si je pridobil bogato teoretično in praktično znanje, ki ga je sedaj pripravljen deliti tudi z nami. Obema novima članoma Uredniškega odbora se zahvaljujemo za njuno sodelovanje pri reviji Acta Geotechnica Slovenica. Bojana Dolinar Glavna urednica EDITORIAL The Editorial Board of Acta Geotechnica Slovenica is pleased to announce that the number of articles received by the journal has increased significantly. This demonstrates the improved recognition of our journal and a greater appreciation among our authors. We are, of course, very satisfied with such a response, which enables us to select some really excellent articles. Unfortunately, however, this means we have to reject a lot of high-quality papers due to the current limitations of our journal. Therefore, only about ten per cent of all the received papers can be published, and the authors of other quality articles are recommended to look to other journals. Last year we already increased the number of published articles from ten to twelve, and in 2015 we have published thirteen papers. Consequently, the publication costs are higher, which is a big problem for the University of Maribor, Faculty of Civil Engineering, Transportation Engineering and Architecture, because co-financing from the Slovenian Research Agency has been drastically decreased. In spite of this, we do not want to pass on the costs to the authors, putting them in an unequal position for publication. Last year we decided to admit new members to our Editorial Board. Above all, the selection criteria consider the quality of previous reviews and the individual's willingness to be involved in future cooperation. Therefore, Dr. Ana Petkovšek from the University of Ljubljana, Faculty of Civil and Geodetic Engineering, Chair of Soil Mechanics with Laboratory joined us last year because of her experiences in geomechanics, road building and environmental geotechnics. Dr. Kosta Urumović from the Croatian Geological Survey, Department of Hydro-geology and Engineering Geology joined us this year. His work with the Croatian Geological Survey mostly focuses on hydrology, the field in which he has gained a great deal of theoretical and practical knowledge, and which he is willing to share with us. We want to thank both new members of our Editorial Board for their cooperation with Acta Geotechnica Slovenica. Bojana Dolinar Editor-in-chief NUMERICNO MODELIRANJE MEHANSKEGA OBNAŠANJA INTAKTNE KAMNINE Z MODELOM KONTINUUMA IN Z MODELOM Z VORONOIJEVIMI ELEMENTI - OBČUTLJIVO-STNA ANALIZA Izvleček У zadnjih letih je numerično modeliranje mehanskega obnašanja intaktne kamnine na mikroskopskem nivoju doseglo veliko pozornosti med številnimi raziskovalci. Eden od možnih načinov predstavitve mikrostrukture intaktne kamnine so modeli s poligonalnimi elementi, ki imajo ploskovne kontakte. У tej raziskavi se je mikrostruk-turo intaktne kamnine generiralo v programskem okolju UDEC (Universal Distinct Element Code) z uporabo Voronoijevega mozaika s poligonalnimi elementi. Njihovo mehansko obnašanje definiramo z mikroparametri, ki jih ne moremo neposredno ugotoviti s standardnimi laboratorijskimi preiskavami, temveč moramo izvesti numerično kalibracijo intaktnega materiala. Postopek kalibracije poteka tako, da iščemo ustrezen niz mikroparametrov modela, ki čim bolje opiše mehansko obnašanje kamnine v laboratoriju. S tem namenom smo izdelali numerične modele za Brazilsko preiskavo, direktno natezno preiskavo, enoosno tlačno preiskavo in dvoosno preiskavo. Z izvajanjem simulacij v omenjenih modelih je bila narejena obsežna parametrična občutljivostna analiza, ki je podala boljše razumevanje vpliva vhodnih mikroparame-trov na odziv modela ter je predvidela relacije med mikro-parametri in makroskopskim odzivom modela. Rezultate teh analiz lahko uporabimo kot splošne smernice za bolj učinkovito in kvalitetno izvedbo kalibracije intaktnega materiala. Paralelno z Voronoijevim modelom smo izvajali še numerične simulacije s kontinuum modelom, kjer je bilo konstitutivno obnašanje modela opisano z Mohr--Coulombovim modelom. Izkazalo se je, da Voronoijev model bolje napove razmerje med Brazilsko in direktno natezno trdnostjo, nastanek razpok in način porušitve ter rezidualno obnašanje modeliranega materiala. Teja Fabjan (vodilni avtor) prej pri IRGO Consulting, d.o.o. Slovenčeva 93, 1000 Ljubljana, Slovenija E-pošta: tejafabjan@gmail.com Diego Mas Ivars Itasca Consultants AB Isafjordsgatan 39B, SE-164 40 Kista, Švedska E-pošta: diego@itasca.se Vladimir Vukadin Inštitut za rudarstvo, geotehnologijo in okolje Slovenčeva 93, 1000 Ljubljana, Slovenija E-pošta: vlado.vukadin@irgo.si Ključne besede metoda diskretnih elementov, parametrična občutljivostna analiza, intaktna kamnina,Voronoijev mozaik, mikromehanski parametri, standardne laboratorijske preiskave NUMERICAL SIMULATION OF INTACT ROCK BEHAVIOUR VIA THE CONTINUUM AND VORONOI TESSELLATION MODELS: A SENSITIVITY ANALYSIS Keywords distinct-element method, parametric sensitivity analysis, intact rock, Voronoi tessellation, micromechanical properties, standard laboratory test Abstract The numerical simulation of intact rock microstructure and its influence on macro-scale behaviour has received a lot of attention in the research community in recent years. Generating a grain-like structure with polygonal area contacts is one of the avenues explored for describing the rock's microstructure. A Voronoi tessellation implemented in the Universal Distinct-Element Code (UDEC) is used to generate models with a polygonal microstructure that represent intact rock. The mechanical behaviour of the Voronoi polygons is defined by micro-properties, which cannot be measured directly in the laboratory. A numerical calibration procedure is needed to produce the macroscopic response of a model that corresponds to the material behaviour measured during a laboratory experiment. In this research, Brazilian, direct tensile, uniaxial compressive and biaxial test models are constructed to simulate the intact rock behaviour under a standard laboratory stress. An extensive series of parametric sensitivity analyses are executed in order to understand the influence of the input micro-properties on every model test behaviour and predict the relation between the micro-properties and the model's macro response. The results can be treated as general guidelines for a complete and efficient intact rock calibration procedure. In parallel, a continuum-based model using the Mohr-Coulomb constitutive relationship is running as a benchmark. It has been shown that the Voronoi-based models through their microstructure approach better reproduce the Brazilian to direct tensile strength ratio, and show a better representation of the dilation, crack pattern and post-peak behaviour in comparison to continuum models. Teja Fabjan (corresponding author) Formerly at IRGO Consulting, d.o.o. Slovenčeva 93, 1000 Ljubljana, Slovenia E-mail: tejafabjan@gmail.com Diego Mas Ivars Itasca Consultants AB Isafjordsgatan 39B, SE-164 40 Kista, Sweden E-mail: diego@itasca.se Vladimir Vukadin Institute for Mining, Geotechnology and Environment Slovenčeva 93, 1000 Ljubljana, Slovenia E-mail: vlado.vukadin@irgo.si 1 INTRODUCTION It is well known that even apparently intact rock is heterogeneous on the scale of the grains [1-3]. Grain shape, grain minerals and different types of other defects included in the matrix introduce a heterogeneity in the internal stress distribution [4]. This has a direct influence on the mechanical behaviour of the intact material in terms of fracture initiation and fracture growth, which control the rock's failure, deformability and crack pattern under loading [5, 6]. All these features control the mechanical behaviour of intact rock. In the past two decades numerical simulations of intact rock microstructure have received a lot of attention from the research community. Many different models have been developed, ranging from the early models based on the continuum finite-element method, to the more recent distinct- or particle-element method models. Continuum-based models (e.g., RFPA, [7]; R-T2D, [8]; DIGS, [4]) can simulate the rock failure process by introducing the heterogeneity of rock properties and using the principles of linear elastic fracture mechanics (LEFMs). The medium is assumed to be composed of many mesoscopic elements with different material properties [9]. Local mechanical parameters, both in terms of elastic and strength properties, are reduced following a stochastic distribution, e.g., a Weibull distribution [3, 9]. However, the choice of input properties can be subjective and highly dependent on the statistical distribution of the parameters [3]; furthermore, the complete detachment of the square-shaped elements is truncated. Since Potyondy & Cundall [10] showed that rock-like material can be successfully modelled with the Bonded Particle Model (BPM), a significant improvement was achieved in modelling the microstructure of rock material. The BPM uses the distinct-element method, as implemented in PFC [11], where a continuum constitutive model is replaced by a simple particle/contact logic [12, 13]. The BPM assumes that the rock behaves like a cemented granular material, where rigid circular or spherical particles are bonded together at their contact points [10]. Reaching the desired material macro-scale behaviour is achieved through a calibration of the input micro-properties [14, 15, 16]. However, when the material is calibrated to the uniaxial compressive strength, the standard BPM significantly over predicts the tensile strength, gives very low triaxial strengths and the failure envelope is linear [17, 18]. By introducing a clumped-particle geometry, or particle cluster logic, in combination with the BPM, these limitations are significantly reduced [10, 18]. Even with the use of clumped-particles and clusters the BPM still suffers from a limitation to match the low ratio of tensile to compressive strength [6], the lowest reachable being at 1/14 [18]. Ratios reported in the literature have much larger variations, ranging from 1/50 to 1/3, depending on rock type, porosity, grain size and other heterogeneities in the rock microstructure [19, 20]. Significant contributions were made by introducing the Grain-Based Model (GBM), where circular/spherical-shaped particles are grouped into polygons with angular boundaries [6]. In this manner, polygonal grains resemble the mineral structure as observed in crystalline rock [12]. Nevertheless, GBM can still not match the Brazilian strength-to-compressive-strength ratio observed in the laboratory [21]. Recently, the BPM has been further developed by incorporating the flat-joint model for grain contacts, allowing for the introduction of a polygonal grain structure that improves the tensile-to-compressive-strength ratio [22]. On the other hand, polygon-shaped grains can be generated directly by using the Voronoi or Dalaunay tessellation generators already implemented in the Universal Distinct-Element Code (UDEC) [23]. Polygonal random-shaped particles are interlocked with area contacts that are defined by well-known geomechanical input properties (e.g., cohesion, tensile strength, friction angle, etc.). Lan et al. [2], Kazerani & Zhao [12, 24], Alzo'ubi [25], Kazerani [13] and Gao [26] reported promising results in terms of reaching a reasonable tensile-to-compressive-strength ratio, curved failure envelope, common failure pattern observed in rocks, etc. The purpose of our study was to evaluate the potential of the Voronoi tessellation model for the simulation of intact rock behaviour. Therefore, a detailed parametric study of the input parameters and their influence on the model's macro behaviour is carried out. The mechanical behaviour of a generic intact rock is simulated by the distinct-element method using UDEC. The rock material is subjected to four standard laboratory test models: Brazilian, direct tensile, uniaxial compressive and biaxial. In parallel, a continuum-based model using a Mohr-Coulomb constitutive relation is running as a benchmark. 2 NUMERICAL MODELLING Numerical distinct-element method simulations are performed using the two-dimensional software UDEC [23]. In UDEC, the mechanical behaviour of a discontinuous system is represented by (a) the behaviour of discontinuities (contacts) and (b) the behaviour of solid material (grains, polygonal particles, blocks). The contact behaviour is described by any joint model, while blocks may be assumed as rigid or deformable. If deformable blocks are chosen, consisting of triangular finite-strain elements (zones), a constitutive behaviour should be defined by using any available constitutive continuum model or user-defined model. As with any other discrete element code, it allows finite displacements and rotations of discrete blocks, including complete detachment, and recognizes new contacts automatically as the calculation progresses [23]. 2.1 Time-marching algorithm UDEC follows an explicit time-marching scheme that solves the equation of motion directly and is identical to that used by the explicit finite-difference method for continuum analyses. The movement disturbance of a discontinuous system is caused by the applied loads or body forces. At the start of each time step the procedure first recognizes the new contact locations from known, fixed, block positions. The force-displacement law is then applied to each contact that generates new contact forces. Based on the resultant force and the moment at the block centroids arising from the known forces acting on them, Newton's second law is applied to all the blocks. The velocity and angular velocity are updated and the linear and angular accelerations are calculated. At this point a calculation step could be repeated by calculating the new block location and the block rotation based on the known velocity and angular velocity. If the blocks are deformable, the motion is calculated at the grid points of triangular elements within the block, and the block material constitutive relations give new stresses within the elements [11, 23]. One time step is needed for every cycle around the loop. The time step is limited by the assumption that velocities and accelerations are constant within the time step. To avoid dynamic effects and to reach quasi-static conditions, each time-step increment should be sufficiently small so that disturbances cannot propagate between one discrete element and its intermediate neighbours. The procedure is repeated until either a satisfactory state of equilibrium or one continuing failure results [23]. 2.2 Constitutive behaviour of intact rock UDEC allows for many possibilities to represent block behaviour by using any of the already-implemented constitutive models. In this study the intact rock is modelled with both continuum and Voronoi tessellation approaches. In the continuum approach the intact block behaviour is assumed to follow the conventional Mohr-Coulomb constitutive model, and is from here on denoted as the continuum model. In summary, the general concept follows a set of linear equations in the principal stress space, a1, a2, a3, where the shear failure of the isotropic elastic material is controlled by the cohesion c and the friction angle ф in the shear yield function: fs = s1 -< 1 + sinj 1 — sinj + 2cj1+^ (1) 1 — sinj Tensile failure is controlled by the tension yield function, where the minor principal stress a3 cannot exceed the tensile strength at: s t ft = s3 — a (2) More details about the Mohr-Coulomb failure criterion can be found in the software manuals [23] or elsewhere (e.g., [27]). On the other hand, the rock microstructure is modelled as an assemblage of distinct deformable polygons, and is subsequently denoted as the Voronoi model. UDEC has an in-built, automatic generator of the Voronoi tessellation pattern, where a particular region in a model can be subdivided into randomly sized polygons. In such a manner, polygons can represent grains, assemblies of grains and/or other flaws (i.e., defects) in the intact rock that are randomly disturbed over the sample [4, 23] (Fig. 1). a) b) c) d) Figure 1. Schematic presentation of Voronoi tessellation generator logic used in UDEC: (a) random generation of points controlled by seed number, (b) generation of Delaunay triangulation, (c) generation of Voronoi tessellation and (d) Voronoi polygons model block [28]. The Voronoi algorithm begins by distributing points within the tessellation region (Fig. 1a). The seed number should be inserted as an input value to set the seed for the random-number generation used by the Voronoi tessellation. During this step the points are allowed to be moved according to the iteration procedure. A larger number of iterations will generate a more uniform spacing between the points and therefore a more uniform tessellation will be generated. These points are then interconnected with lines that form triangles, called Delaunay triangulation, and thus each triangle is circumscribed within a circle containing its three vertexes (Fig. 1b). Finally, the Voronoi polygons are created by constructing perpendicular bisections of all the triangles that share a common side (Fig. 1c). The polygons are truncated at the boundaries of the tessellation region (Fig. 1d) [23]. The Voronoi polygons are assigned an isotropic elastic deformable material model, which means that failure on the micro-scale only takes place at the polygon boundaries, not inside them. The Voronoi contacts are assumed to behave following the Coulomb slip area joint model. The stress-displacement relation in the normal and shear directions is assumed to be linear and governed by the normal stiffness kn and the shear stiffness ks such that: Aan =-knAUn (3) Ats = -ksAues (4) where Aan is the effective normal stress increment, Ars is the effective shear stress increment, Ли/ is the elastic component of the incremental shear displacement and ЛиГ1 is the normal displacement increment. However, the normal stress is limited by a limiting tensile strength (if an < -at , then an = 0) and the shear stress is limited by a linear combination of the cohesion c and the friction angle ф: Ы £ c + sn tanj = Tmax (5) When the shear stress exceeds the maximum shear stress (| Ts I > Tmax) then: T s = sign(Aus )rmax (6) where Лus is the total incremental shear displacement. In addition, joint dilation may occur at the onset of the slip of the joint. The accumulated dilation f is limited either by a high normal stress level or a large accumulated shear displacement that exceeds a limiting value umax. Therefore, if |Ts| < Tmax or I Ts I = Tmax and I us I > umax, then f = 0, otherwise f > 0. 3 NUMERICAL LABORATORY The characterisation of intact rock is achieved by performing standard laboratory tests, so the aim of this study is to recreate the same testing procedure within a numerical environment using the distinct-element method. Therefore, four standard laboratory test models are developed: a splitting tensile strength test (Brazilian test), a direct tensile test, a uniaxial compressive strength (UCS) test and a biaxial test. These tests characterise the rock material in the tensile and compressive stress paths and can be used to calibrate the constitutive behaviour of the modelled material to the real material. 3.1 Model generation Two types of model geometry are generated. A circle-shaped model is used to simulate the Brazilian test and a rectangular-shaped model is used to simulate the direct tensile test, the UCS test and the biaxial test. Plane stress analyses are adopted for all the cases to assume thin model samples. These models are run by applying a vertical velocity at the external boundary at the top of the sample. The velocity is fixed and set to zero in the vertical direction at the external boundary at the bottom of the sample. The only difference between the tests is the direction of the applied velocity. The velocity is directed upwards in the direct tensile test to produce tensile conditions in the sample, while in the other test models the velocity is directed downwards to generate a compressive stress state. In addition, to ensure confinement in the biaxial test model, a horizontal compressive stress is acting on the left and right external boundaries, and the vertical horizontal compressive stress acts on the top external boundary of the model sample. The boundary conditions are shown in Fig. 2. Figure 2. Boundary conditions in (a) Brazilian test, (b) direct tensile test, (c) uniaxial compressive test and (d) biaxial test. "F" denotes fixed grid points and "S" denotes confinement stress boundary conditions. Numbers on samples show locations of monitoring grid points. 3.2 Monitoring The monitoring procedure is very similar in all the test models, where the grid-point forces, the displacement and the positions of the grid points are monitored using a FISH function. Based on these monitored data points the strength and deformation properties are calculated for every time step while running the test analysis. The applied vertical velocity at the top of the sample causes a reaction force at every grid point in the sample. In each time step the total sum of the vertical forces Fy for every grid point at the external boundary on the top side of the sample is divided by the cross-section d of the sample to obtain the vertical stress oy. The peak stress is determined as the failure strength of the tested material. However, in the Brazilian test model the reaction force is converted into stress by using the standard Brazilian strength ffß equation [29]: = 2Fy/7 dt (7) where Fy is the total sum of the vertical forces, d is the sample (circle) diameter and t is the sample thickness, which is assumed to be unity in a 2D analysis. Material deformability is calculated by monitoring the displacements and the positions at the grid points in the middle of every side of the sample, as shown in Fig. 2. The horizontal deformation is calculated based on monitoring locations 1 and 3, and the vertical deformation is calculated based on monitoring locations 2 and 4. The Young's modulus is then calculated by dividing the axial stress by the vertical deformation, and the Poisson's ratio is calculated by dividing the horizontal deformation by the vertical deformation. 4 NUMERICAL PARAMETRIC SENSITIVITY ANALYSIS_ An extensive series of parametric sensitivity analyses are executed in order to understand the influence of the input properties on the model's behaviour and the relationship between the micro-properties and the model's macro response. This kind of study has given us a set of general guidelines for the calibration procedure of the intact material. Parallel to this, a series of analyses on continuum models are carried out to compare the results with those of the Voronoi model. The parametric sensitivity analysis is conducted by running a series of simulations in which each significant input parameter is varied systematically over a chosen range, while all the other parameters were fixed. Some of the ranges used here are not realistic for rock-like material. However, as will be shown later, some valuable insights into the numerical environment can be obtained from them. Two types of parameters were tested for all four numerical laboratory tests: (a) numerical input properties (particle size, seed number, number of iterations, mesh density, number of fixed grid points, strain rate, model sample size and shape) and (b) material input properties (material density, Young's modulus, Poisson's ratio, normal and shear stiffness, cohesion, friction angle, tensile strength and dilation). The input properties for both the continuum and the Voronoi models are arbitrarily chosen within a range that can represent the behaviour of sedimentary rock, such as marl, siltstone, etc. The input properties are listed in Table 1 and Table 2 (next page) for the continuum and Voronoi models, respectively. Table 1. Input properties used in the parametric sensitivity analysis for the continuum model. Parameter Input value of parameter Tested range of parameter Numerical input properties Model height (mm) 500 50-1,500 Model width/diameter (mm) 250 25-750 Zone edge length (mm) 10 10-4,000 Applied velocity (m/s) 0.001a 0.00-0.3 Material input properties Material density (kg/m3) 2,548 10-100,000 Young's modulus (GPa) 8.940 1-100 Poisson's ratio ( ) 0.167 0.1-0.4 Cohesion (MPa) 2.2 0.5-20 Friction angle (°) 30 0-80 Limit tensile strength (MPa) 1.5 0.1-20 Dilation angle (°) 2 0-15 a Applied velocity for UCS test and biaxial test was 0.005 m/s. 4.1 Size of the Voronoi polygons The resolution of the model is controlled by the size of the Voronoi polygons and it affects the peak strength and the material deformability. The size of the Voronoi polygons is defined by the Voronoi edge length, which must be at least 20 times that of the corner rounding length and specifies the input value of the average edge length of the Voronoi polygons [23]. A higher model resolution requires more cycles (time) to reach mechanical equilibrium, and also needs more time to generate the initial polygon packing in a model than a lower model resolution. Thus, a series of analyses are run to determine the optimal model resolution. A minimum of ten seed distributions within any test model were analysed. As will be explained later, different seed numbers give Table 2. Input micro-properties used in the parametric sensitivity analysis for the Voronoi model. Parameter Input value of parameter Tested range of parameter Numerical input properties Model height (mm) 500 50-25,000 Model width/diameter (mm) 250 25-12,500 Voronoi edge length 16.5 6.25-80 Zone edge length (mm) 8.25 3-125 No. of iterations ( ) 80 13-1000 Seed no. ( ) 222 10-999 Applied velocity (m/s) 0.001b 0.001-0.01 Material input properties Material density (kg/m3) 2,548 10-100,000 Young's modulus of polygons (GPa) 9.2 1-100 Poisson's ratio of polygons ( ) 0.2 0.1-0.4 Normal stiffness (GPa/m) 10,000 100-100,000 Shear stiffness (GPa/m) 2,000 20-200,000 Contact cohesion (MPa) 2.2 0-10 Contact friction angle (°) 10 0-80 Contact limit tensile strength (MPa) 1.5 0.1-20 Dilation angle (°) 2 0-15 b Applied velocity for UCS test and biaxial test was 0.005 m/s. different peak strengths, so the mean peak strength and the relative standard deviation (RSD) are calculated from the results of the tests. In order to give a better representation of the results from all the model tests, the peak strength values are normalized to the mean value. Figure 3. Effect of polygon size resolution on normalized peak strength and corresponding RSD error for the Voronoi model. The results for all four tests shown in Fig. 3 imply that at least 10 Voronoi polygons along the cross-section of a model are required to minimize the polygon resolution's influence on the peak strength to an acceptable level. In this range the RSD begins to stabilize at around 5 %. However, both the uniaxial compressive and biaxial model tests show a slight increase in the peak strength when the number of Voronoi polygons is larger. That effect is marginal and might be a result of statistical repetition. When more Voronoi elements are included, the fracture (i.e., the failure surface) can have a different path and the resulting peak may be slightly different. Based on this, the recommended model resolution is between 10 and 25 Voronoi polygons per model sample cross-section. A denser Voronoi size resolution can be used, but it will lead to increased calculation times. According to the ASTM standard [30] and ISRM suggested methods [31], the sample diameter tested in the laboratory should be larger than the largest mineral grain by at least 10 or 20 times. If it is assumed that the Voronoi polygons represent grains, this will be in agreement with the presented results. However, the results from the Brazilian strength (Fig. 3) do not show the same trend as the other model tests. The strength keeps decreasing, even though a maximum of 40 Voronoi polygons per model test cross-section are tested. The reason for the decreasing trend of the Brazilian strength might be due to the fact that there are fixed grid points at the external boundary, which are located on a curved line due to the circular model geometry. The same response is observed in the Brazilian model tests with the Dalanuay triangles [12]. However, to keep the consistency in the results from every model test in the numerical laboratory, the polygon resolution was kept equal to 15 from this point on in the study for all the samples tested. The same as in the BPM models [10], the particle size resolution slightly influences the elastic properties. Increasing the polygon resolution up to 40, the Young's modulus decreases and the Poisson's ratio increases, both by approximately 15 % (Fig. 4). A possible explanation for this behaviour might be as follows. While the polygon resolution is increasing, more contacts are formed in the model. Therefore, the model sample has a greater ability to deform due to the displacement along the contacts, which cause a lower Young's modulus and a higher Poisson's ratio. 4.2 Microstructure heterogeneity The heterogeneity of intact rock at the grain scale is due to the grain shape and size, the grain contact length and orientations, the constituent minerals and other 9600 - LEGEND » Direct tensile test -*- Uniaxial compressive test » Biaxial test 7200 - 5 10 15 20 25 30 35 40 a) No. of polygon;; in cosrì. scdion i : 45 0.23 ? 0.22 - 0.21 0.2 ■ 0.19 LEGEND Direct tensile test ♦ Uniaxial compressive test * Biaxial test 5 10 15 20 25 30 35 b) No. of polygons i- cross sco'.io" f ; Figure 4. Effect of polygon size resolution on (a) Young's modulus and (b) Poisson's ratio for the Voronoi model. spatial defects [2]. In the Voronoi model the material's geometric heterogeneity and the polygon distribution can be controlled by the seed number and the number of steps in the iteration procedure applied in the Voronoi tessellation generation. Every seed number generates a random array of Voronoi tessellation and thus every tested model will give a slightly different macro response, e.g., every seed number will produce a different peak strength, crack pattern, etc., as can be seen from Fig. 5. As previously mentioned, using the recommended polygon resolution (10 to 25) an approximately 5 % RSD error in the peak strength is expected by testing the model with a different seed number. Similarly, the variability of the results due to microstructure heterogeneity occurs in the laboratory as well, because in reality every sample has a slightly different microstructure. Therefore, the distributions of results can be compared to the variations observed in the laboratory tests. However, the major difference between the numerical and laboratory results is in the size of the variability of the results. In the laboratory the size of the variations is greatly affected by the grain size and other micro-cracks in the rock. This means that some rocks will have a smaller variation than others [28, 32]. On the other hand, as can be seen from Fig 3., the RSD error (which is basically the size of the variations of the results) is going to stabilize when the number of Voronoi polygons is more than 10 polygons per sample cross-section. Therefore, in numerical simulations, following the approach proposed in this paper, we will get approximately the same RSD error, no matter what kind of rock we would like to model. On the other hand, the number of iterations controls the Voronoi tessellation uniformity, as can be seen from Fig. 6. Increasing the number of iterations reduces the heterogeneity, modifies the crack pattern and influences the peak strength. In general, increasing the iteration Figure 6. Example of uniaxial compressive test model sample at failure tested with the same seed number, but different iteration number (a) 13, (b) 150 and (c) 800. Figure 7. Effect of number of iteration steps on normalized peak strength for the Voronoi model. number increases the peak strength in all the model tests. The peak strength tends to stabilize using a large number of iteration steps. In tensile tests (e.g., the direct tensile and Brazilian tests) stabilisation is reached after 150 iterations, but in compressive tests (e.g., the UCS and biaxial test) stabilization is not reached until 550 iteration steps are applied (Fig. 7). This response can be directly linked to the tessellation geometry, where the difference in the geometry becomes more and more insignificant, when comparing two Voronoi tessellations generated with a very high iteration number (higher than 600 iteration steps). Moreover, increasing the number of iteration steps has a negligible influence on the peak-strength scatter between the different seed numbers, and over the entire tested range this is approximately 5 % of RSD. In most cases, rocks have a heterogeneous microstructure geometry [2]. This can be controlled with the number of iterations. According to this study, less than 250 iterations are recommended to produce natural material behaviour. Iteration numbers between 50 and 250 generate a relatively homogeneous model microstructure, which can be compared to sedimentary rocks (Fig. 6b). A small number of iterations (below 50) generates a more heterogeneous model microstructure, similar to magmatic rocks, e.g., granites (Fig. 6a), but gives a lower peak strength. Finally, a very high iteration number (more than 250) can produce unnatural crack patterns that are normally not observed in real rock (Fig. 6c). 4.3 Mesh density The blocks/polygons in both the continuum and Voronoi models are discretized into deformable, triangular, finite-difference zones and its zone edge length is one of the input parameters for the model (Fig. 8). To achieve reasonable calculation times and stable numerical results, an optimal zone edge length (mesh density) must be defined for the present model size. For that reason a sensitivity analysis of the zone edge length's influence on all the model tests is carried out. Figure 8. Schematic representation of Voronoi polygon blocks discretized into deformable triangular finite-difference zones. The sensitivity analysis with the continuum model has shown that the influence of mesh density within the range considered in this study of the peak strength and deformability can be neglected in rectangular-shaped samples (Fig. 9), e.g., direct tensile test, uniaxial compressive test and biaxial test, but it is evident in circular model samples, e.g., the Brazilian test model (Fig. 10). Since the mesh density does not affect the peak strength and deformability in the rectangular model samples within the range considered in this study, the same mesh density can be used, regardless of the model size. However, increasing the model size with a fixed zone edge will affect the peak strength in the Brazilian test model. Therefore, the results for all the test models will be consistent if the appropriate mesh density for the Brazilian test model is found and then used in other test models. Figure 9. Effect of zone edge length on normalized peak strength for the continuum model. Figure 10. Dependence between Brazilian strength and model sample size to zone edge length ratio for four model sizes (diameter) for the continuum model. As shown in Fig. 10, stabilization of the Brazilian strength occurs in the range of model size to zone edge length ratio between 25 and 35 (with a scatter error in the range of 5 %). Higher ratios lead to a longer calculation time, with only small improvements in the model prediction. The reason for the decreasing trend of the Brazilian strength was already discussed in Section 4.1. On the other hand, the mesh density slightly influences the material behaviour in the Voronoi model, as shown in Fig. 11. The normalized strength is dependent on the ratio between the Voronoi edge length and the zone edge length. Since the zone edge length is defined as the maximum edge length of the triangular zones, there will be a critical limit with regard to the Voronoi edge length (i.e., the zone edge length is smaller, or equal, to the Voronoi edge length). When the ratio is less than one, there is no influence on the peak strength, regardless of the increasing value of the zone edge length, because the zone edge length is now limited by the Voronoi edge length. Figure 11. Effect of Voronoi edge length to zone edge length ratio on normalized peak strength for the Voronoi model. With the ratio increasing from 1 to 4, the peak strength decreases by about 5 % for all the tests, except for the direct tensile test, where it increases by the same amount. A possible reason for this could be that the direct tensile test has a tensile stress state that is in contrast to the other three tests, where we have predominantly compressive stress states. Also, the peak strength scatter between the different seed numbers at an arbitrary mesh density is approximately 5 % of RSD. Based on these findings and to avoid long calculation times, the appropriate ratio of the Voronoi edge length to the zone edge length is between 1.5 and 4. 4.4 Number of fixed grid points As discussed in Section 3.1, the velocity is applied to grid points (called fixed grid points) on the top side of the sample. In the rectangular-shaped model samples the number of fixed grid points is related to the mesh density in the continuum model and to both the mesh density and the polygon size in the Voronoi model. Therefore, the number of fixed grid points is a consequence of the mesh and polygon density and has no direct influence on the model's behaviour. On the other hand, the circular geometry of the Brazilian test model contains several edges along the outer boundary. Besides the mesh density and polygon size, the number of edge segments is related to the number of fixed grid points in the upper and lower external boundaries. Increasing the number of segments implies an increase in the number of fixed grid points in both models. The results of the sensitivity analysis show that an increase in the number of fixed grid points will increase the Brazilian strength until stabilization is reached (Fig. 12 on next page). An up to 10 % RSD error can be expected in the stable region. Therefore, at least 4 fixed grid points (Fig. 12b) or 80 segments (Fig. 12a) are recommended values for both the continuum and Voronoi models. f g 'S SE (D 0.9 0.8 0.7 0.6 0.5 0.4 03 0.2 LEGEND ■*■ Continuum model ■*■ Vororioi model 50 100 150 200 250 a) V; . 300 £ 1.1 - 1 0.9 o.e 0.7 0.6 0.5 0.4 03 0.2 А LEGEND Continuum model Voronoi model 1 2 3 4 5 6 7 8 b) IV ;.: ■ l.;;';!i - ; : 10 Figure 12. Effect of (a) number of segments and (b) number of fixed grid points on the Brazilian strength for the Voronoi model. 4.5 Strain rate In accordance with the ISRM [33], a quasi-static laboratory compression test is achieved within the strain-rate interval 10-5 to 10-4 /s. Faster or slower strain rates than these can have a significant impact on the material's behaviour [34]. A similar response is observed in numerical simulations. Running a model test with a high applied velocity introduces dynamic effects in the simulation. However, choosing a very low strain rate is inefficient and leads to very time-consuming calculations. The aim of this part of the sensitivity analysis is to find the optimal strain rate for the current modelling study. The optimal applied velocity can be determined using a trial-and-error procedure based on a comparison between the different stress-strain curves obtained from tests running at different velocities. As the maximum applied velocity for which the shape of the stress-strain curve stabilizes and the peak strength becomes constant, the applied velocity can be adopted as the optimal velocity for the current model analysis. However, dealing with velocity is not very convenient when trying to compare the results between the different model tests. In that case it is better to use the displacement rate Лер which can be calculated based on the number of cycles Ncyc needed for the current model analysis at the chosen applied velocity Dil and the mechanical time step ktm taken by UDEC: De = Dil, ■btm/N, cyc (8) Several Voronoi model sizes with the same mesh density and model heterogeneity are tested. The results from the uniaxial compressive test simulations presented in Fig. 13a show that there is a relationship between the applied velocity in the Voronoi model and the size of the model. LEGEND Velocity 0.01 m/s (Straight) ■ a Velocity 0.005 m/s (Straight) - ft ■ Velocity 0.001 m/s (Strenght) * - Vekraty 0.01 m/s (Displ.rate) Velocity 0.005 m/s (Displ.rate) Velocity 0.001 m/s (Displ.rate) 0.1 1 10 a) 1E-003 1.2 1.15 1.1 1.05 1 0,95 -0-90.350.80.75 LEGEND Brazilian test Direct tensile test Uniaxial compressive test Biaxial test 1E-006 1E-005 0.0001 b) Figure 13. (a) Effect of applied velocity on normalized peak strength and displacement rate for different model sizes in uniaxial compressive test in the Voronoi model, and (b) effect of applied velocities on normalized peak strength in the continuum model. Larger model sizes require a lower applied velocity than the smaller model sizes. That model behaviour is a result of the numerical algorithm of the distinct element code. Therefore, from that graph (Fig. 13a) an optimal displacement rate can be determined, regardless of the model size. For example, the peak strength stabilization of the solid blue line (an applied velocity of 0.01 m/s) is reached at a model size of approximately 0.5 m. Dragging a vertical line from here onto the dashed blue line shows that stabilization happens at approximately 10-5 mm/step, which can be adopted as an optimal displacement rate. A similar procedure can be used for an applied velocity of 0.005 m/s (red line). The peak strength of the tested sample with an applied velocity of 0.001 m/s is stable over the entire tested interval, because the displacement rate is lower than the optimal displacement rate. To conclude, the optimal displacement rate should be 10-5 mm/step for all the model tests in this study. A similar displacement rate can be adopted for the continuum model (Fig. 13b) and for all the model tests performed in the numerical laboratory (Section 3). 12000 11000 10000 - 0.01 a) _oy - Applied ve:ouily in/s: 0.23 0.22 0.15- 0.14 LEGEND Direct tensile test Uniaxial compressive test Biaxial test 0.001 b) Loc 0.01 ■ Applied velocity (m/s) The applied velocity also has a small influence on the Young's modulus and the Poisson's ratio. As can be seen in Fig. 14: at a high applied velocity the Young's modulus is higher, but the Poisson's ratio is lower. This can be interpreted as follows. At a high applied velocity (displacement rate) the model becomes stiffer (a higher Young's modulus), because it has less ability to accommodate all the deformation (a lower Poisson's ratio). On the other hand, at a low displacement rate, when the time step is small enough for the model to accommodate more of the deformation (higher Poisson's ratio), it becomes softer (lower Young's modulus). With a further decrease in the displacement rate the Young's modulus and the Poisson's ratio tend to stabilize. 4.6 Model shape aspect ratio It is well documented that the sample shape influences the material strength and the ductility. As the aspect ratio of the height to width (h:w) increases, the strength and ductility decrease by decreasing the horizontal confinement stress in unconfined tests [35]. Due to the friction between the sample ends and the platen, and the difference between the elastic properties of the rock and the platen, the sample will be restrained near its ends and prevented from deforming uniformly [34]. Different rectangular-shaped model aspect ratios are tested under direct tensile, uniaxial compressive and biaxial test conditions. The results from the Voronoi model show that the peak strength decreases with an increasing aspect ratio (Fig. 15). In addition, the post-peak response becomes more ductile with a decreasing aspect ratio (Fig. 16). The peak strength stabilizes with an aspect ratio higher than 2, which is in accordance with the standard recommendations (i.e., a suggested Figure 14. Effect of applied velocity on (a) Young's modulus and (b) Poisson's ratio for the Voronoi model. Figure 15. Effect of model shape aspect ratio (h:w) on the normalized peak strength for the Voronoi model. Lateral / Axial strain ( ) Youna's modulus (input) (GPa) Figure 16. Uniaxial compressive test model stress-strain curves Figure 17. Effect of input Young's modulus on the calculated for different values of height to width (h:w) aspect ratio for the Poisson's ratio for the Voronoi model. Voronoi model. aspect ratio between 2.0 to 2.5 [36], and 2.5 to 3.0 [33]). As expected, the model aspect ratio has a negligible effect in direct tensile tests. 4.7 Material density According to UDEC theory the internal density is irrelevant to the modelling in static systems since the internal forces in the system are low [23]. The material density only influences the calculation time, because it is proportionally related to the critical time step [23]. We have confirmed this for both the continuum and Voronoi models [28]. Figure 18. Effect of Young's modulus on normalized peak strength for the Voronoi model. 4.8 Deformability properties: Young's modulus, Poisson's ratio, normal and shear stiffness The input Young's modulus and Poisson's ratio have a direct influence on the calculated Young's modulus and the Poisson's ratio in the continuum and Voronoi models. However, the input Young's modulus also influences the calculated Poisson's ratio in the Voronoi model as well. As can be seen in Fig. 17, an approximately 40 % change can be expected in the tested range. The explanation may be that a higher input Young's modulus causes the polygons to behave in a stiffer manner, so the majority of the deformation is accumulated by the contacts and thus the material's Poisson's ratio increases. When the input Young's modulus is lower, the polygons are softer and can accumulate the majority of the deformation, and thus the value of the calculated Poisson's ratio becomes lower. It is also shown that the Young's modulus has a negligible influence on the peak strength in all the model tests, except for the Voronoi model Brazilian test, where an increasing Young's modulus decreases the Brazilian strength by about 40 % in the tested range (Fig. 18). In the Voronoi model the deformability is controlled by two additional micro-properties, i.e., the contact normal stiffness kn and the contact shear stiffness ks. However, both have an influence on the strength behaviour of the material as well. The effect of the normal and the shear stiffness are shown (Fig. 19) separately in terms of the normalised strength and stiffness ratios kn/ks and ks/kn. Different shapes of the peak-strength trend curves can be distinguished between the peak strengths measured in the tensile tests (e.g., the direct tension test and the Brazilian test) and the compressive tests (e.g., the uniaxial compressive test and the biaxial test). The compressive strength is slightly increasing throughout the tested interval, while the tensile strength quickly stabilizes. Observations of 1.4 1.3- € 1.2 1.1 - 1 - 0.9 0.8 10 kn ( ) 20 30 40 50 ks / kn О LEGEND -*- Brazilian test - Direct tensile test Uniaxial compressive test ■ Biaxial test —è- 20 40 60 80 a) --■:> iii::l -t ""■■ .Л ■:: - : 100 1.6 1.4 Dl 1,2 С E 0.8 о 0.6 0.4 12 16 20 .i.i.i. ^_____ ______1 LEGEND -*- Brazilian test ■*■ Direct tensile test *■ Uniaxial compressive test ■*■ Biaxial test 40 80 120 b) i . i'i. 160 200 Figure 19. Correlation between normalized peak strength to (a) normal stiffness and kn/ks ratio, and to (b) shear stiffness and ks/kn ratio for the Voronoi model. the model's material behaviour (e.g., the shape of the stress-strain curve) is in agreement with the intact rock behaviour known from the laboratory if the stiffness ratio kn/ks is approximately between 1 to 20 (i.e., the stiffness ratio ks/kn is approximately between 1.0 and 0.05). The empirical equation developed by D. Martin and reported by Vallejos et al. [14] and the data reported by Li & Wong [37] show that the laboratory-measured Brazilian strength is consistently higher than the direct tensile strength [28], most likely due to the friction between the platen and the sample ends, which causes confinement stress in the sample. The aim of our study is not to quantify the strength ratio, but only to check the trends observed in the laboratory, i.e., whether the simulated Brazilian strength is lower, higher or the same compared to the direct tensile strength. It is found that the laboratory-measured strength ratios between the Brazilian tests and the direct tensile tests can be matched by choosing an appropriate stiffness ratio kn/ks. It is difficult to define a reasonable interval, because it is dependent on other input micro-properties as well, but is in between kn/ks = 1 to 20 range. In fact, the tensile-to-compressive strength ratio is also influenced by the stiffness ratio kn/ks and by other micro-properties (e.g., cohesion, friction angle, and tensile strength limit). The normal to shear stiffness ratio can control the Young's modulus and the Poisson's ratio, as has been already recognized by Kazerani & Zhao [24] and in the BPM model [10, 17]. When the stiffness ratio is higher than a critical value, the calculated elastic properties become constant, even if the normal and/or the shear stiffness changes (Fig. 20, Fig. 21). Therefore, the stiff- 12.000 6.000 -4.000 2.000 0 10 kn /ks( ) 20 30 40 50 ______ -- LEGEND Direct tensile test ■+■ Uniaxial compressive lesi Biaxial test 20 40 60 80 a) 100 ks / kn () a 12 1B 20 11,000 - a, ooo ■ é 7,000 E P 6,000 5,000 4,000 3,000 H ----1 LEGEND » Direct tensile test ■*■ Uniaxial compressive test * Biaxial test i . i . i 40 80 120 160 b) 200 Figure 20. Correlation between Young's modulus to (a) normal stiffness and kn/ks ratio, and to (b) shear stiffness and ks/kn ratio for the Voronoi model. 10 кг I ks ( ) 20 30 40 50 r LEGEND Direct tensile test Uniaxial compressive test Biaxial test 20 40 60 00 a) I'm-"-i bT'in^s .'1 li r>." 100 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 ks I kn ( ) 8 12 16 20 LEGEND Direct tensile test Uniaxial compressive test Biadai test 40 80 120 160 b) si n-^-sfj ■.'С '^mj 200 Figure 21. Correlation between Poisson's ratio to (a) normal stiffness and kn/ks ratio, and to (b) shear stiffness and ks/kn ratio for the Voronoi model. ness ratio can be calibrated to some extent to fit the material's Poisson's ratio and the Young's modulus. 4.9 Strength properties: cohesion, friction angle and tensile strength limit In Voronoi models, micro-properties, such as the contact cohesion, contact friction angle and contact tensile strength limit of the Voronoi polygonal contacts, exert the main control over the material's strength behaviour. In the continuum model the material's strength behaviour is controlled by the (macro) cohesion, the (macro) friction angle and the (macro) tensile strength limit. This is confirmed by the series of sensitivity analyses presented in this section. First, the sensitivity analysis on cohesion is discussed. As shown in Fig. 22 a similar model response can be observed in both the continuum and the Voronoi models. The uniaxial or biaxial compressive strength show a linear dependence on increasing cohesion, which is controlled by the shear yield function. However, the Brazilian and direct tensile strength are controlled either by the tension or the shear yield function, as shown in Fig. 23. At first, when the cohesion is lower than the tensile strength limit (which can be an unrealistic case), it shows a linear dependence on the tensile strength and is controlled by the shear yield function (Fig. 23a). Once the tensile strength limit is reached (which is a more realistic case), the calculated peak tensile strength becomes constant, regardless of the cohesion increase and is controlled by the tension yield function (Fig. 23c). LEGEND Brazilian test Direct tensile test Uniaxial compressive tesi Biaxial test 6 б 10 12 14 16 18 20 a) LEGEND Brazilian test Direct tensile test Uniaxial compressive test Biaxial test 23456769 10 b) Figure 22. Correlation between cohesion and normalized peak strength for (a) the continuum and (b) the Voronoi model. Oj~ (TT > G GT > Gg G G aL7" G Figure 23. Schematic presentation of Mohr circles at (a) (absolute) high tensile strength limit, (b) similar values of cohesion and tensile strength limit, (c) high cohesion and (d) high friction angle. Moreover, a relatively high cohesion compared to the tensile strength limit could produce undesirable post-peak behaviour in terms of a strong material resistance during sliding, which generates residual stresses that are commonly higher than the peak strength (mostly in the Brazilian test model). In this sense, the material model's behaviour is in agreement with the intact rock behaviour known from the laboratory if the cohesion values are approximately up to five times the tensile strength limit. Similar to the sensitivity analysis of cohesion, the continuum and the Voronoi model show comparable results with respect to the friction angle (Fig. 24). The uniaxial and biaxial compressive strength show an exponential dependence on the increasing friction angle, which is controlled by the shear yield function. The Brazilian and direct tensile strength, on the other hand, is controlled either by the tension or the shear yield function. At low and moderate friction angles the direct tensile strength is constant and is controlled by the tensile yield function (Fig. 23c). However, a small decrease in the Brazilian strength is observed at low friction angles, which is more obvious in the continuum model (Fig. 24a). At present, such a response is not fully understood. On the other hand, very high friction angles (which may not be a realistic case) induce a decrease in the tensile strength as well, but this time the shear yield function controls the response (Fig. 23d). The sensitivity analysis has also shown that high friction angles can induce a strong shear resistance in the post peak, so that the residual shear strength exceeds the peak shear strength (Fig. 25). If such a post-peak behaviour is undesirable, it can be limited by choosing lower friction angles as the input (micro) property. The major difference between both models is seen in the parametric sensitivity analysis of the tensile strength limit. The continuum model shows a very similar response to that of the sensitivity analysis of cohesion, except that the (uniaxial or biaxial) compressive strength is unaffected by the tensile strength limit (Fig. 26a). The Brazilian strength and the direct tensile strength exhibit a clear dependence upon the tensile strength limit (Fig. 26a). At first, when the tensile strength limit is low, the (macro) tensile strength is LEGEND Brazilian test Direct tensile test Uniaxial compressive test Biaxial tesi 20 30 40 50 a) 20 30 40 50 60 b) Figure 24. Correlation between friction angle and normalized peak strength for (a) the continuum and (b) the Voronoi model. /I / I , _ / / I /''-""'' / i „--z*-:...,-- LEGEND — S 0' 40° ■ 80° 0 0.001 0.002 a) 0.003 m 0.8 to 0. S 06 j= ■s 0,4 ì 0.2 0.001 0.002 0.003 0.0W 0,005 b) Figure 25. Effect of friction angle on stress-strain curves of Brazilian test model for (a) the continuum and (b) the Voronoi model. LEGEND Brazilian Lest Direct tensile test Uniaxial compressive test Biaxial test 2 4 6 Э 10 12 14 16 a) 18 20 2 4 6 Э 10 12 14 16 b) Figure 26. Correlation between tensile strength limit and normalized peak strength for (a) the continuum and (b) the Voronoi model. controlled by the tension yield function (Fig. 23c). After reaching the critical value, e.g., the maximum tensile strength (otmax = c/tanp), the tensile strength becomes constant and the failure is controlled by the shear yield function (Fig. 23a). On the other hand, the Voronoi model shows a semilinear dependence on the tensile strength limit (Fig. 26b). Increasing the magnitude of the tensile strength limit increases the failure strength in every model tested. At this point it is clear that the macro-strength behaviour of the Voronoi model is not directly controlled by the input value of each single micro-property (i.e., it needs to be calibrated as a whole because all the micro-properties interact in complex ways) unlike the continuum model. In fact, every micro-property has an influence on the model's response, and has to be taken into account during the material's calibration procedure. Contrary to what is observed in the Voronoi model, the Brazilian strength can never overcome the direct tensile strength in the continuum model. The difference between the tensile strengths in the continuum model is mainly affected by the cohesion and the tensile strength limit. The direct tensile strength is higher than the Brazilian strength as long as the Mohr-Coulomb model does not include the tensile strength limit. This is satisfied when the tensile strength limit is much higher than the cohesion (which is basically an unrealistic case, see Fig. 23a) or when the cohesion and tensile strength limits are similar in magnitude (Fig. 23b). When the cohesion is significantly larger than the tensile strength limit (Fig. 23c), the Brazilian strength and the direct tensile strength are similar in magnitude. (It should be noted that the graphs in the paper do not show the absolute, but the normalized, strength values.) The Brazilian strength could slightly overcome the direct tensile strength in some simulations due to the rough mesh density (Section 4.3) or an inappropriate number of circle segments (Section 4.4). The Brazilian strength to direct tensile strength ratio in the Voronoi model is described in Section 4.8. 4.10 Dilation Since intact rock is represented as an assembly of polygons with micro strength properties at their contacts, the cracking and failure that occurs during loading can potentially reproduce macro dilation and control the material's macro strength. For the Voronoi model, a sensitivity analysis on input dilation angle, which is assigned to contacts, shows no significant influence on the pre-peak region and the peak strength in both models (Fig. 27). After failure, when the polygons start to slide along their contacts, this causes the material to dilate so that a higher dilation angle produces a lower residual strength (Fig. 27b). As expected, dilation has no effect on direct tensile-test simulations [28], since polygons are only moving apart. On the contrary, in the continuum model, a higher dilation angle produces a higher residual strength (Fig. 27a). The reason is in the formulation of the Mohr-Coulomb constitutive model, where material with a larger dilation angle tends to follow a strain-hardening path. 5 CONCLUDING REMARKS The numerical modelling of intact rock microstructures has received a lot of attention in the past two decades. It has been recognized that rock microstructure is of key importance in intact rock behaviour. For this purpose a Voronoi polygonal-shaped model was used to generate the intact rock microstructure and simulate its response to four standard laboratory tests. A detailed parametric sensitivity analysis was performed to obtain an in-depth understanding of the model behaviour. The findings from the parametric sensitivity analysis can be summarized as follows: - The effect on the peak strength of the Voronoi polygon size resolution can be eliminated by having more than 10 polygons in the model sample's cross-section. - Voronoi polygon heterogeneity and distribution can be controlled by the seed number and the number of iterations. Up to 250 iteration steps are necessary to produce a natural crack pattern, which is also different for different seed numbers. - Voronoi polygons are discretized into deformable zones. The recommended ratio of Voronoi edge length to zone edge length is at least 1.5. - The strain rate effect is affected by the model size and in this project it is eliminated when it is lower than 10-5 mm/step. - The model shape height-to-width ratio affects the peak strength and the material ductility, which is in accordance with the laboratory experiments. - The stiffness ratio kn/ks affects the rock peak strength, the Young's modulus and the Poisson's ratio. Reasonable material behaviour can be achieved if the stiffness ratio kn/ks is between 1 and 20. The stiffness ratio kn/ks control affects the Brazilian to direct tensile strength ratio as well. 0.0003 0.0006 0.0009 a) 8 7 6 £L c sl to й л щ 4 У) ™ 3 X < 2 1 0 0.001 0.002 0.003 b) 0.004 Figure 27. Effect of dilation on the stress-strain curves of the uniaxial compressive test for (a) the continuum and (b) the Voronoi model. Contact cohesion, contact friction angle and contact tensile strength limit affect the material's tensile and compressive peak strengths. Contact dilation affects the post-peak material beha- In parallel with the Voronoi model, the continuum model was running as a benchmark, as it is widely known and used in many engineering applications. Since the Voronoi model can better represent the intact rock microstructure it was recognized that it could be a better approach to study damage processes in rocks during laboratory testing. Being able to control the model microstructure is one of its advantages. By using different seed numbers in the Voronoi models' generation, a slightly different crack pattern for the same conditions is reached. This makes it possible to match the model's macro response to real material behaviour not only quantitatively, but also qualitatively. This is not possible with the continuum model, because of its limitations as an elasto-plastic model and its linear failure envelope. Moreover, Voronoi models better reproduce the Brazilian to direct tensile strength ratio, and show a better representation of the dilation and post-peak behaviour in comparison to the continuum models. Acknowledgements The operation part of this research is financed by the European Union, European Social Fund. The authors also wish to thank the Itasca Education Partnership Program, Itasca Internationl, Inc. and the Geotechnical Laboratory at the Institute for Mining, Geotechnology and Environment (IRGO, Slovenia). Mr. Abel Sanchez Juncal from Itasca Consultants AB is also greatly acknowledged for his contribution to various parts of this study. REFERENCES_ [1] Marin, D. 1993. The strength of massive Lac du Bonnet granite around underground openings. PhD thesis. Canada: University of Manitoba. [2] Lan, H., Martin, C.D., Hu, B. 2010. Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading. J. Geophys. Res. 115(B1):B01202. [3] Mahabadi, O.K., Randall, N.X., Zong, Z., Grasselli, G. 2012. A novel approach for micro-scale characterization and modeling of geomaterials incorporating actual material heterogeneity. Geophys. Res. Lett. 39, L01303. [4] Van de Steen, B., Vervoort, A., Napier, J.A.L., Durrheim, R.J. 2003. Implementation of a flaw model to the fracturing around a vertical shaft. Rock Mech. Rock Eng. 36, 143-161. [5] Chen, S., Yue, Z.Q., Tham, L.G. 2007. Digital image based approach for three-dimensional mechanical analysis of heterogeneous rocks. Rock Mech. Rock Eng. 40, 145-168. [6] Potyondy, D.O. 2010. A Grain-Based Model for Rock: Approaching the True Microstructure. Proceedings of Rock Mechanics in the Nordic Countries, Kongsberg, Norway, pp. 10. [7] Tang, C.A. 1997. Numerical simulation on progressive failure leading to collapse and associated seis-micity. Int. J. Rock Mech. Min. Sci. 34(2), 249-261. [8] Lui, H. 2003. Numerical Modelling of the Rock Fracture Process under Mechanical Loading. PhD thesis. Sweden: Lulea University of Technology. [9] Wang, S.Y., Sloan, S.W., Huang, M.L., Tang, C.A. 2011. Numerical Study of Failure Mechanism of Serial and Parallel Rock Pillars. Rock Mech. Rock Eng. 44, 179-198. [10] Potyondy, D.O., Cundall, P.A. 2004. A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci. 41, 1329-1364. [11] Itasca Consulting Group, Inc. 2008. PFC 3D, Particle Flow Code in 3 Dimensions. 4th ed. Minneapolis, USA. [12] Kazerani, T., Zhao, J. 2013. A Microstructure-Based Model to Characterize micromechanical Parameters Controlling Compressive and Tensile Failure in Crystallized Rock. Rock Mech. Rock Eng. 47(2), 435-452. [13] Kazerani, T. 2013. A discontinuum-based model to simulate compressive and tensile failure in sedimentary rock. J. Rock Mech. Geotech. Eng. 5, 378-388. [14] Vallejos, J.A., Brzovic, A., Lopez, C., Bouzeran, L., Mas Ivars, D. 2013. Application of the Synthetic Rock Mass approach to characterize rock mass behavior at the El Teniente Mine, Chile. Continuum and Distinct Element Numerical Modeling in Geomechanics. Minneapolis, Paper 07-02. pp.15. [15] Brzovic, A., Schachter, P., de los Santos, C., Valle-jos, J., Mas Ivars. D. 2014. Characterization and Synthetic Simulations to Determine Rock Mass Behaviour at the El Teniente Mine, Chile. Part I. Proceedings, Caving 2014, 3rd International Symposium on Block and Sublevel Caving, Santiago, Chile, pp. 8. [16] Vallejos, J., Suzuki, K., Brzovic, A., Mas Ivars, D. 2014. Characterization and synthetic simulations to determine rock mass behaviour at the El Teniente mine, Chile. Part II. Proceedings, Caving 2014, 3rd International Symposium on Block and Sublevel Caving, Santiago, Chile, pp. 9. [17] Diederichs, M.S. 1999. Instability of hard rock masses: the role of tensile damage and relaxation. PhD thesis. Canada: University of Waterloo. [18] Cho, N., Martin, C.D., Sego, D.C. 2007. A clumped particle model for rock. Int. J. Rock Mech. Min. Sci. 44, 997-1010. [19] Sheorey, P.R. 1997. Empirical rock failure criteria. Rotterdam: Balkema. [20] Vutukuri, V.S., Lama, R.D., Saluja, S.S. 1974. Handbook on mechanical properties of rocks, testing techniques and results, Vol I. Trans Tech Publications. [21] Bahrani, N., Potyondy, D., Pierce, M. 2012. Simulation of Brazilian Test Using PFC2D Grain-Based Model. Proceedings of 21st Canadian Rock Mechanics Symposium (CARMA), Quebec, Canada, pp. 485-493. [22] Potyondy, D.O. 2012. A flat-jointed bonded-particle material for hard rock. In: Proceedings of 46th U.S. Rock Mechanics/Geomechanics Symposium, Chicago, USA, pp. 10. [23] Itasca Consulting Group, Inc. 2014. UDEC, Universal Distinct Element Code. 4th ed. Minneapolis, USA. [24] Kazerani, T., Zhao, J. 2010. Micromechanical parameters in bonded particle method for modelling of brittle material failure. Int. J. Numer. Anal. Meth. Geomech. 34, 1877-1895. [25] Alzo'ubi, A.K. 2012. Modeling of rocks under direct shear loading by using discrete element method. AHU. J. of Engineering & Applied Sciences, 4(2), 5-20. [26] Gao, F. 2013.Simulation of Failure Mechanisms around Underground Coal Mine Openings Using Discrete Element Modelling. PhD thesis. Canada: Simon Fraser University. [27] Labuz, J.F., Zang, A. 2012. Mohr-Coulomb Failure Criterion. Rock Mech. Rock Eng. 45, 975-979. [28] Fabjan, T. 2015. Numerical modelling of mechanical properties of discontinuities in jointed and heterogeneous rock mass. PhD thesis. Slovenia: University of Ljubljana. [29] Bieniawski, T.Z. 1977. Suggested Methods for Determining Tensile Strength of Rock materials. Int. J. Rock Mech. Min. Sci. Geomech Abs. 15(3), 99-103. [30] ASTM International 2013. Standard test method for compressive strength and elastic moduli of intact rock core specimens under varying states of stress and temperatures. ASTM D7012-13. [31] Fairhurst, C., Hudson, J. 1999. Draft ISRM suggested method for the complete stress-strain curve for intact rock in uniaxial compression. Int. J. Rock Mech. Min. Sci. 36(2), 279-289. [32] Yoshinaka, R., Osada, M., Park, H., Sasaki, T., Sasaki, K. 2008. Practical determination of mechanical design parameters of intact rock considering scale effect. Eng. Geol. 96, 173-186. [33] Bieniawski, T.Z., Franklin, J.A., Bernede, M.J., Duffaut, P., Rummel, F., Horibe, T., Broch, E., Rodrigues, E., Van Heerden, W.L., Vokler, U.W., Hansagi, I., Szlavin, J., Brady, B.T., Deere, D.U., Hawkes, I., Milanovic, D. 1979. Suggested methods for determining the uniaxial compressive strength and deformability of rock materials. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 16(2), 135-140. [34] Brady, B.H.G, Brown, E.T. 2006. Rock Mechanics for underground mining. 3rd ed. Dordrecht: Springer. [35] Hudson, J.A., Harrison, J.P. 2005. Engineering rock mechanics an introduction to the principles. Pergamon; Elsevier Science. [36] ASTM International. Standard Practice for Preparing Rock Core as Cylindrical Test Specimens and Verifying Conformance to Dimensional and Shape Tolerances. ASTM D 4543-08. [37] Li, D., Wong, L.N.Y. 2012. The Brazilian Disc Test for Rock Mechanics Applications: Review and New Insights. Rock Mech. Rock Eng. DOI 10.1007/ s00603-012-0257-7. OCENA STABILNOSTI OBALNIH PEČIN Z UPORABO DIGITALNIH POSNETKOV Izvleček Preiskovano območje okoli naselja Stara Baška (otok Krk, SV del kanala Jadranskega morja) je v občutljivem geodinamičnih ravnovesju. Morska erozija je zelo vidna in med 1966 in 2004 je evidentiran umik pečin med 4 in 5 metrov. Pobočja pečin oblikujejo strme breče. Močni valovi in formacije zarez vsled rezanja valov predstavlja glavni vzrok nestabilnosti. Sekundarni vzroki so preperevanja in erozije zemljine in hribine na pečini. Udor pečine pobočja se lahko pojavi v skalnih gmotah višjih parametrov trdnosti, kjer so zareze pomaknjene nekaj metrov v vznožje pečine pobočja. Preizkušena je bila kombinirana metoda za analizo stabilnosti obalnih pečin, ki vključuje model konzolnega nosilca in SfMfotogrametrijo. S to metodo lahko dobimo zelo detajlirane 3D geometrijske podatke pečine, ki se lahko kasneje uporabijo v izračunih stabil-nostnega modela. To je še posebej pomembno pri stabil-nostnih analizah litološko heterogenih kamninah, kot so breče s spremenljivo geometrijo, ki jih je težko nadomestiti s pravokotnimi površinami. Z enostavno in uporabno SfM metodo presežemo meje tradicionalnega merjenja pri ocenjevanju previsa pečine in dolžine zareze. Igor Ružić University of Rijeka, Faculty of Civil Engineering Radmile Matejčić 3, HR-51000 Reka, Hrvaška E-pošta: igor.ruzic@uniri.hr Čedomir Benac University of Rijeka, Faculty of Civil Engineering Radmile Matejčić 3, HR-51000 Reka, Hrvaška Ivan Marović University of Rijeka, Faculty of Civil Engineering Radmile Matejčić 3, HR-51000 Reka, Hrvaška Suzana Ilić University of Lancaster, Lancaster Environment Centre LA1 4YQ Lancaster, Velika Britanija Ključne besede morska erozija, stabilnost klifov, SfM fotogrametrija, otok Krk, Jadransko morje A STABILITY ASSESSMENT OF COASTAL CLIFFS USING DIGITAL IMAGERY Keywords marine erosion, cliff stability, SfM photogrammetry, Krk island, Adriatic sea Abstract The investigated area around the Stara Baška settlement (the island of Krk, NE channel part of the Adriatic Sea) is in a delicate geodynamic balance. Marine erosion is quite prominent and the recorded cliff retreat between 1966 and 2004 was from 4 to 5 metres. The cliff slopes are formed in talus breccias. Strong waves and formations of wave-cut notches are the main causes of the cliff's instability. The secondary causes are the weathering and erosion of the soil and rocks on the cliff. The slump of the cliff slope can occur in a rock mass with higher strength parameters, where the notches are cut a few metres inward into the toe of the cliff slope. A combined method for the stability analysis of the coastal cliffs was tested; this incorporates the cantilever-beam model and the structure-from-motion (SfM) photogrammetry. This method can provide highly detailed 3D geometrical data of the cliff, which can then be used in the calculations of the stability model. This is particularly important in a stability analysis of lithologically heterogeneous rocks such as breccias with varied geometry, which cannot be easily replaced by a rectangular surface. The simple and useful SfM method overcomes the limitations of traditional surveys in estimating the cliff overhang surface and the notch length. 1 INTRODUCTION The occurrence and extent of coastal cliff erosion are primarily influenced by the lithology, the rock mass structure, and its geotechnical properties. These factors affect the local rates of retreat and the resulting cliff morphology (Grigss and Trenhaile, 1997). These have been the subjects of a number of studies. However, the erosion of hard rock cliffs has received less attention in previous studies due to their relative stability and generally slower rates of retreat than those found in soft rock cliffs. The recession of hard rock cliffs is usually the result of wave undercutting and the subsequent removal of debris deposited at the toe of the cliff by marine processes (Trenhaile, 2002). The most effective erosional processes take place near the water level, leading to the development of cliff notches at the base of the cliffs. These are then found to be responsible for cliff failure due to the imbalance of the forces and moments in the upper cliff-face material (Davies et al., 1991; Williams et al., 1993; Hungr et al., 2001; Matsukura, 2001; Castedo et al., 2012). Two main types of cliff mass movement are identified: topples and slides (Sunamura, 1992). Sub-aerial processes controlled by local climates such as seepage and surface run-off erosion can intensify the decay of the active coastal cliffs and hence their instability. The cantilever-beam model has been used for a stability analysis of coastal cliffs with the presence of notches in a Igor Ružić University of Rijeka, Faculty of Civil Engineering Radmile Matejčić 3, HR-51000 Rijeka, Croatia E-mail: igor.ruzic@uniri.hr Čedomir Benac University of Rijeka, Faculty of Civil Engineering Radmile Matejčić 3, HR-51000 Rijeka, Croatia Ivan Marović University of Rijeka, Faculty of Civil Engineering Radmile Matejčić 3, HR-51000 Rijeka, Croatia Suzana Ilić University of Lancaster, Lancaster Environment Centre LA1 4YQ Lancaster, United Kingdom number of studies (e.g., Thorne and Tovey, 1981; Abam, 1997, Matsukura, 1988, 2001; Kogure et al., 2006; Kogure and Matsukura, 2010; Castedo et al., 2012). The model is based on the assumption that the cliff-overhanging material above the notch acts as cantilever-beam load. The tensile stresses, which are distributed on the upper section, and the compressive stress, distributed on the lower section of the cliff face, are derived from the model (Timoshenko and Gere, 1972). To calculate the stresses, input parameters such as the rock density and the geometry of the cliff face are required. So far the model was mainly used to estimate the instabilities of the slopes in cliffs formed in jointed limestones (Kogure, 2006). In that case the cliff and the fallen block cross-sections have a roughly rectangular shape, which can be simply represented by a simple rectangular geometry (Kogure and Matsukura, 2010). The cliff's cross-section is calculated by multiplying the average cliff height by the horizontal distance between the retreat point of the notch and the averaged location of the cliff face. In the case where the cliffs were formed in lithologically heterogeneous rocks, such as breccias, the cliff profiles have an irregular shape and no visible mechanical discontinuities and hence the geometry parameters, such as averaged cliff-face location, are difficult to define. Some authors (Williams et al., 1993; Kogure and Matsukura, 2010; Castedo et al., 2012) argue that a precise coastal 3D topography is needed for the cliff-stability analysis and the safety-factor determination. Field measurements of the 3D topography are extremely difficult to obtain and in some cases almost impossible due to the restricted access to the site and the risk of a cliff collapse. Moreover, the notch geometry is often difficult to measure with traditional instruments due to its substantial lengths and low heights. Also, a GPS-RTK device often does not work in such confined spaces. Remote sensing methods that overcome these problems have been mainly focused on analysing historical aerial photos and satellite images until recently. These, together with a traditional survey of the location of cliffs tops, have been used to estimate the amount and the cliff-retreat rate (Kuhn and Prüfer, 2014). However, these methods are not sufficient to derive the cliff geometry necessary for the stability models described above. The new methods for topographic-data acquisition, such as LIDAR and those based on structure-for-motion (SfM) photogrammetry (James et al., 2013), enable the collection of high-accuracy and high-density 3D spatial data for researching the coastal geomorphol-ogy. The advantage of the SfM method in comparison with LIDAR, and other traditional photogrammetric methods, is that the images can be taken with a readily available, single digital camera such as 10-megapixel Ricoh GR Digital IV. Multiple photographs are taken from various horizontal and vertical angles (James and Robson, 2012). Unlike traditional photogrammetric methods, which need the 3D location and position of the camera(s), the SfM method solves the problems of camera position and scene geometry simultaneously and automatically by using a highly redundant bundle adjustment based on matching features in multiple overlapping, offset images (Snavely et al., 2008; Westoby et al., 2012). It has been demonstrated that this simple method is highly applicable for studies of gravel beach dynamics (James et al., 2013) at the Uboka pocket beach (Fig. 1, location 2), changes in beach slopes at the gravel beach and erosional processes in talus breccia (Ružić et al., 2013) in Mošćenička Draga (Fig. 1, location 3) and cliff geomorphology (Ružić et al., 2014) at the site of Stara Baška (Fig. 1, location 1). Figure 1. Study area with the location of Stara Baška (1); locations of previous researches in Uboka (2), and Mošćenička Draga (3) and the location of mareograph (tidal gauge) Bakar (4). The north-eastern Adriatic coast is predominantly a steep, rocky coast formed of karstified carbonate rocks, compared to the flat south-western Italian side. For this reason, the geomorphology of this coastal area is more complex and diverse, comprising a number of coastal forms on a small scale (Pikelj and Juračić, 2013). The coastal cliffs near Stara Baška in the Kvarner Bay area in Croatia (Fig. 1, location 1, Fig. 2) are made of talus breccias and have pronounced notches. These cliffs are very unstable and have significantly retreated during the past five decades and collapsed rocky blocks are often found on the beaches fronting the cliffs. Some of the island's impor- tant infrastructure, such as the tourist camp in Stara Baška, is threatened by these cliff collapses. Hence, a prediction of this type of geohazard is extremely valuable for future spatial planning and development on this coastline. The aim of this paper is to examine the stability of the breccias cliffs at the site in Stara Baška using the cantilever-beam model (Kogure et al., 2006) in combination with the SfM method, which is used to derive the geometrical properties of the cliffs. We examine the effect of notches on the distribution of the stresses inside the cliffs as well as on the spatial distribution of the stresses along the cliff face. The advantage and limitations of both methods are illustrated as well as the recommendation for further studies. 2 BACKGROUND TO THE INVESTIGATED SITE The area around the Stara Baška settlement has a delicate geodynamic balance, and the coastal cliff's erosion is intense. The coastal geomorphology, shown in Figure 2, consists of several headlands and embayments. The figure also shows the settlement of the village itself and the camping site on the top of the cliffs. 2.1 Geomorphology and geological setting The Kvarner area, a semi-enclosed section of the Adriatic Sea, is located between the Istrian peninsula and the Vinodol-Velebit coast. Geographically, the island chains Cres-Lošinj and Krk-Rab-Pag divide it into the Rijeka Bay; the Kvarner Bay, Kvarnerić; and the Velebit Channel (Fig. 1). Intensive morphogenetic processes caused by tectonic movements and rapid sea-level changes, as well as climatic changes, resulted in the present shape of the Kvarner area, as well as the shape of Krk Island (Benac and Juračić, 1998; Benac et al., 2013). In the terrestrial part of the Kvarner area, sedimentary carbonate rocks prevail, whereas siliciclastic rocks are rare. Pleistocene and Holocene deposits locally cover this bedrock substrate. The waves in the channel part of the northern Adriatic Sea are not as high as in the western open zone due to the relatively short wind fetch. On a carbonate rocky coast bio-erosional processes prevail over mechanical erosion and tidal notches are common features of these processes. In places where the rock mass is tectonically crushed or karstified, wave-notches and cliffs evolved. Marine erosion and more rapid cliff recession are obvious in a less-resistant siliciclastic rocky coast and Pliocene-Quaternary sediments (Juračić et al., 2009). The area surrounding Stara Baška consists of Upper Cretaceous and Palaeogene carbonate rocks. Palaeogene siliciclastic rocks such as marls and flysch outcrop can be found in a relatively narrow coastal zone. They are mostly covered by Pliocene-Quaternary sediments (Benac et al., 2013). The carbonate rocks, such as limestone, dolomitic limestone and dolomitic breccia, have a different degree of fissuring and form a typical bare karstic landscape. Research was focused on two embay-ments carved into Pliocene-Quaternary talus breccias (Fig. 2). The outcrops of flysch rocks are partially visible in the narrow coastal zone. Site research showed that breccias have pronounced horizontal stratification and the joints were not visible. Horizontal layers consist of limestone fragments/clasts from a few millimetres in size to blocks greater than 50 Figure 2. Simplified engineering geological map of the investigated area (red rectangle): 1-carbonate rocks (Upper Cretaceous and Palaeogene limestones), 2-siliciclastic rock mass (Palaeogene marls and flysch) mostly covered by Quaternary talus breccias, 3-Quater-nary alluvial sediments, 4-periodical surface water flow, 5-partially canalized water flow. cm, with prevailing fragments of 1-4 cm in size. The matrix is reddish silt-sandy cement (20-50% clay) with a different degree of calcification. A different degree of calcification has a great influence on the erodibility as well as on the strength parameters of breccias. Cliffs formed in talus breccias have an average height of 5 to 6 metres. The top layer at the subsurface, of about 1 m thickness, is completely weathered and has soil features. These features are clearly visible in the block fallen from the cliff, as shown in the photo taken during one of the field visits (Fig. 3). On the beach in front of the cliffs, gravel and partially coarse sand sediments prevail. The average gravel grain size is between 2 and 20 mm. north-eastern bora forms moderate waves due to small fetches, despite having the highest speed. South-east winds generate the highest waves, as shown in Figure 4. Figure 4 shows significant wave-height simulations using the SWAN model (Booij et al., 1999). The wave simulation was performed for a SE wind speed of 25 m/s (100-year return period). For illustration purposes, the model assumes a uniform distribution of wind over the domain and does not take into account the complex wind field in the Bay. The area around the Stara Baška settlement has an approximately 25% lower significant wave height than the exposed opposite channel part (Fig. 4). Figure 3. Recently collapsed cliff (09/01/2014), site location 1 (SL1): 1 - soil (like weathered zone), 2 - talus breccia. Figure 4. Significant wave height (in metres), NE part of Kvarnerić (according to the SWAN wave model). 2.2 Climatic conditions Low and moderate winds intersected with silent periods are most frequent, while storm winds (speed > 30 m/s) are rare (Leder et al., 1998). North-eastern wind bora is predominant and reaches a maximum speed of 30 m/s. The mean sea level recorded at the nearest tidal gauge in Bakar (Fig. 1, location 4) is 0.15 m (vertical datum of the Republic of Croatia - CVD). This is a micro-tidal environment with a tidal range between 30 and 35 cm (Benac et al., 2004). The Kvarner area is a semi-enclosed section of the Adriatic Sea between the Istrian peninsula and the Vinodol-Velebit coast. This affects the increase in water levels and wave heights during storms. The storm surge can increase the water level up to 1.18 m CVD, as recorded in Bakar on 1 December 2008. The 100-year return and the 1-year return periods of the storm surge levels are 1.30 m CVD and 0.80 m CVD, respectively (Ružić, 2003). The wave heights in the Kvarnerić area near Stara Baška are smaller than in the western open zone due to the relatively short wind fetch. Hence, the 2.3 Evidence of cliff retreat The analysis of historical aerial images showed that the breccias are prone to rapid cliff recession (Ružić et al., 2014). The changes in the cliff-top position were identified by a comparison of geo-referenced aerial images from 1966 and ortho-rectified photos from 2004, shown in Figure 5. Autodesk RasterDesign software was used for the georefer-encing. The georeferenced map's estimated horizontal accuracy is up to 0.50 metres due to the distortion of the old areal images. Visible changes were noticed between these two aerial surveys. Furthermore, large rocky blocks situated at the toe of the cliffs indicate earlier local block falls. In the western embayment, a significant cliff retreat of up to 5.5 m is found only in the eastern part (Fig. 5c, near profile 2). This resulted in an increase of the beach area by about 75 m2. The cliff line in the middle part of the embayment (Fig. 5c, profile 1) did not change significantly. In this case, the detected changes are smaller than the possible errors in the image analysis. In the Figure 5. Aerial photographs showing considerable coastal changes between 1966 and 2004: a) geo-referenced aerial photo from 1966; b) ortho-photo map from 2004; c) coastal changes between 1966 and 2004: 1 - coastline 2004; 2 - coastline 1966; 3 - cliff line 2004; 4 - cliff line 1966. Red line and numbers (Fig. 5c) indicate analysed profiles. eastern embayment, the retreat in the cliff-top location is up to 5 meters at the western side (Fig. 5c, profiles 4 and 5). The cliff near profile 4 has retreated by about 2 metres since 1966. The beach area increased here up to 265 m2. This analysis showed that the cliff's retreat results in an extension of the beach area. Also, it showed that there is a noticeable spatial difference in the cliff retreat, which will be explained later. The stability of the cliffs at the locations identified here will be further investigated by using the combined cantilever-beam model and the SfM method. 3 METHODS For the stability analysis, the cantilever-beam model given by Kogure et al. (2006) is adopted and briefly summarized here. Kogure et al. (2006) assessed a cliff-failure process based on three stages: (1) a developing of notch forms as a column like a cantilever beam; (2) this column generates tension cracks on the top surface of the cliff due to the tensile stress arising from its own weight; and (3) a vertical extension of the cracks causes bending failure of this column-shaped mass, which results in a toppling mode of failure. Point-load tests were carried out using a Digital rock strength index apparatus (Control 45-DO55Ü/D). The density in dry and wet conditions of the selected samples in talus breccia has been measured according to standard methods (HRN EN 1097-6). The estimation of the Geological Strength Index (GSI) has provided us with the according standard methods (Hoek and Brown, 1991; Hoek, 1996; Marinos and Hoek, 2000). The maximum bending stress (ffmax) inside a cliff is estimated as the ratio of the bending moment (M) and the modulus of the cliff section (Z) given in Equation 1. M (1) Z 1 2 c"n M = 2 p gbhcl z=1 bh2 6c (2) (3) Equations 2 and 3 are used to estimate the bending moment and the modules for a cliff section with breadth (b), height (hc) and notch depth (ln), respectively. The cliff height (hc) is defined as the vertical distance between the retreat point of a notch (the deepest point of a notch) and the top surface of a cliff. The notch depth ln is defined as the horizontal distance between the retreat point of a notch and an averaged cliff face. For cliffs with tension cracks present on the top of the cliff, only the height that experiences the tensile and comparative stresses would be taken into account. Figure 8 shows the stress distribution inside a cliff with no tension cracks (Fig. 8a) and with tension cracks (Fig. 8b). A structure-from-motion (SfM) photogrammetric method was used to obtain parameters such as the area of the cliff-overhang material, the height of the cliff face and the length of the notch (Ružić et al., 2014). These are derived from geo-referenced 3D point clouds generated from approximately 200 images using the Autodesk ReCap online service (https://recap360.autodesk.com/). The images were collected using a single 10-megapixel Ricoh GR Digital IV camera equipped with a high-contrast 6-mm f/1.9 GR lens, and focusing on embay-ments as shown in Figure 5b (i.e., SL1 and SL2). Special care was needed to acquire sharp and well-focused images from a variety of locations covering a range of vertical and horizontal angles. Also, visual checks were used to select images to upload to the Autodesk ReCap online service and for 3D point cloud verification, as suggested by Podobnikar (2009). Cloud Compare software (http:// www.danielgm.net/cc/) was used for the point cloud geo-referencing. A small number of ground-control points (GCPs) were measured in the field using Real Time Kinematic GPS (RTK-GPS) after the images were taken. These were then used for the transformation of coordinates from a relative to an absolute coordinate system (e.g., Westoby et al., 2012). It is important to point out that all the images were taken within a period of less than two hours. Complete field measurements, including RTK-GPS surveys, took about six hours for a two-person team. 4 RESULTS_ 4.1 Properties of the material in the cliff The stability analysis was based on laboratory-derived data. The material bulk density in dry conditions is 24.5 kN/m3, and 25.3 kN/m3 in wet conditions. The uniaxial comprehensive strengths are USC = 29.21-70.21 MPa in dry conditions and USC = 19.82-62.10 MPa in wet conditions. The estimated Geological Strength Index is GSI = 25-35, and the internal friction angle is F = 28-35° (mi = 20) (Hoek and Brown, 1991; Hoek, 1996; Marinos and Hoek, 2000). It is difficult to estimate the geological parameters in litho-logically very heterogeneous rocks like breccia. A correlation between the heterogeneous breccias compressive and tensile strengths cannot be unambiguously defined. 4.2 Cliff geometry derived from the photogrammetry An example of a high-density 3D point cloud, derived from a set of photo images taken on the sites SL1 and SL2 are shown in Figure 5. Different parts of a cliff and beach, such as the cliff face, notches, vegetation on the top of the cliff, the upper and lower beach, and their structures are clearly identifiable from the cloud. The locations of the cliff and notches and their geometrical properties are obtained from a geo-referenced 3D point cloud. The spatial variation in the cliff geometry is clearly illustrated in Figure 6. The high-resolution of the images and the colour provide additional information about the spatial variation of the cliff geomorphology and the vegetation cover. Cliff cross-sections obtained from a geo-referenced point cloud at a number of locations marked in Figure 6 are shown in Figure 7. All the profiles have vertical convex cliff faces fronted by a gently sloping beach. They can be categorized as those that do not have notches (Profiles 1 and 4) and those with pronounced notches and a cliff overhang (Profiles, 2, 3 and 5). For Profile 1 the cliff base is at 1.8 m above sea level and a half metre above the 100-year storm surge level. The beach width taken from the mean sea level to the toe of the cliff is 11.4 m. For Profile 4 the cliff base is 1.2 m above the sea and the beach width is 5.5 m. The cliff notch is not visible from the derived cliff profile. However, the photorealistic 3D point cloud (Fig. 6) indicates that the cliff undercut is filled with sediments. The cliff profiles 2, 3 and 5 have pronounced notches at the cliff toe. The 100-year surge water level reaches the notches at Profiles 2 and 3, but this is not the case with Profile 5, which is about 40 cm above it. The beach width is 5.6 m and 0.8 m in Profile 2 and Profile 3, respectively. The largest width of the beach is 7.7 m, found at Profile 5. 4.3 Results of stability analysis The cross-section geometry derived from the SfM photogrammetry was used to calculate parameters such as the cliff area, the cliff overhanging material weight and the centre of gravity. These are then used in the cantilever-beam model for an estimation of the cliff's maximum bending stress (ffx-max) using equation 1. The wet density of breccias of 25.3 kN/m3 was taken into account. The estimates for profiles 2, 3 and 5 are given in Table 1 (row 11) and Figure 8a illustrates the distribution of stresses for Profile 3. The stresses are calculated by considering the 'proper' area of the profile rather than a rectangular approximation in the momentum equation (2). These are given in Table 1 in the 15th row. All the profiles' calculations of stresses are corrected for the presence of a soil layer on the top of the cliff. As the soil does not support tensile stresses, it leads to an increase of the tensile and compressive stresses in the underlying rock. These increased stresses were calculated following Kogure et al.'s (2010) approach for cliffs with pronounced cracks. The method reduces the critical height of a cross-section Figure 6. Photorealistic 3D point clouds shown from different angles with marked profiles (red lines)- Figure 7. Cross-sections perpendicular to the coastline, profiles 1-5 (Figs.6 and 10); characteristic sea levels: mean sea level (blue) and 100-year storm surge (magenta). for a tensile crack depth. Here, the critical height of a cross-section was reduced by the thickness of a soil layer (Fig. 8b). From the 3D point cloud shown in Figure 6, Figure 8. Stress distribution inside a cliff (profile 3) with no tension crack (a) and with a tension crack - soil layer (b), and the main variables used in the calculation. the depth of the top soil layer was estimated to be 1 m (as shown in Fig. 8), and the critical height was reduced accordingly (Table 1, row 13). 4.4 Distribution of stresses along the cliff face The detailed cliff geometry derived from the point cloud allowed precise geometrical parameters for the spatial distribution of the stresses along the cliff face. The section around Profile 3 was considered, as it has the widest cliff notch and may be considered as a 3D phenomenon. Two cross-sections on each side from the current Profile 3, with a 2 m offset between them, were chosen (Fig. 9a) for a stability analysis. Figure 9b shows the difference in chosen profiles in terms of the notch length and the cross-sectional area. Table 2 gives a summary of the estimates, parameters and calculated values for all the cross-sections. The largest stresses are found between profiles 3-3 (0.352 MPa/m1) and 3-4 (0.356 MPa/m1). The section stresses decrease approximately 16 times between profiles 3-4 and 3-5, the distance of which is only 2 metres. This clearly correlates with the length of the cliff notches, which is the deepest between profiles 3-2 and 3-4. Table 1. The cliff (beam) model input parameters and the resulting sections stress. Profile Property Symbol Unit 2 3 5 Bulk density Y kN/m3 25.30 25.30 25.30 s ie tr Cliff area A m2 16.20 23.04 33.09 e po Notch depth In m 3.73 5.16 6.56 or p al Centroid distance Ig m 1.56 2.21 2.71 r e n Cliff weight G kN/m1 409.90 582.90 837.10 e a Momentum M kNm/m1 640.70 1290.60 2269.20 Cliff height hc m 5.78 5.69 6.92 l t) Depth of tensile crack ahc m 0.00 0.00 0.00 ts 0 <3 3 h (ac Bearing cross-section height hc-ahc m 5.78 5.69 6.92 Centre of masses Z 3 m 5.57 5.39 7.98 .S Tensile stress a kN/m2/m1 115.10 239.30 284.30 -Я « ^ Depth of tensile crack ahc m 1.00 1.00 1.00 а о ^ cac Bearing cross-section height hc-ahc m 4.78 4.69 5.92 ■s « 8 5? ° 2 Centre of masses Z m3 3.81 3.66 5.84 Tensile stress a kN/m2/m1 168.20 352.40 388.50 Profiles 3-1 and 3-5 are the cliff-boundary profiles, the cliff overhanging material area and the section height ratio is restricted and its stresses are limited (0.035 and 0.023 MPa/m1). The tensile stresses are in places larger than 300 kN/m3, which are approaching the critical stresses for these types of rocks. It is also shown that this combined method can identify much closer locations with maximum stresses. The strength of the tensile stresses and their distribution result in a non-uniform cliff recession, also seen from the analysis of the historical aerial photographs. 5 DISCUSSION_ Assuming that the cliff collapse is solely the result of the notch length and amount of material overhanging the notch, the proposed combined method could be useful for estimating the stresses in the cliff face and its stability. In particular, it is shown that the use of the SfM method can examine a number of cliff profiles at different spatial distances in order to identify a critical profile. So far, the stability analysis was based on a limited number of cliff profiles and was very much dependent on the subjective choice of a critical profile when taking measurements in the field. Another advantage is that we can estimate the length of the notch much more easily. As an illustration, the tensile stresses at Profile 5 are calculated by taking into account the 'assumed rectangular cliff area' and three different values for the notch length: the minimum, maximum and mean value, as suggested by Kogureet et al. (2006). Those stresses calculations are made by using the cliff height hc and the notch length ln following Kogure et al. (2006) and Table 2. Tensile-stress calculation parameters from profile 3-1 to 3-5. Profile 3 Property Symbol Unit 3-1 3-2 3-3 3-4 3-5 Bulk density Y kN/m3 25.30 25.30 25.30 25.30 25.30 s ie tr Cliff area A m2 4.40 20.62 23.04 23.15 3.87 e po Notch depth ln m 2.35 5.21 5.16 5.42 2.46 or p al Centroid distance Ig m 0.73 1.97 2.21 2.20 0.68 r e n Cliff weight G kN/m1 111.30 521.70 582.90 585.70 97.90 e G Momentum M kNm/m1 81.00 1027.70 1290.60 1286.70 66.80 Cliff height hc m 4.75 6.29 5.69 5.66 5.22 li n tn Depth of tensile crack ahc m 0.00 0.00 0.00 0.00 0.00 h (ac Bearing cross-section height hc-ahc m 4.75 6.29 5.69 5.66 5.22 1 -3 £ ä Centre of masses Z 3 m 3.75 6.60 5.39 5.33 4.54 .S Tensile stress o kN/m2/m1 21.6 155.8 239.3 241.3 14.7 k n t)n Depth of tensile crack ahc m 1.00 1.00 1.00 1.00 1.00 У S 3 catc Bearing cross-section height hc-ahc m 3.75 5.29 4.69 4.66 4.22 ■s « 8 5? ° 2 Centre of masses Z m3 2.34 4.67 3.66 3.61 2.97 Tensile stress o kN/m2/m1 34.7 220.2 352.4 356.1 22.5 described in the previous section. The results are given in Table 3 (row 6) and there is a significant difference between these estimates. These are then compared with the calculated tensile stresses, taking the measured area into account. The difference is 29% when the mean notch length is taken into account. If the mean notch length ln is overestimated, differences might be even higher, as shown in Table 3. This again highlights the importance of obtaining accurate cliff-geometry parameters. The SfM method can be used to estimate the depth of the top-soil layer on the basis of the colour difference in the 3D point cloud. However, these estimates are more uncertain as they are based on the information from the cliff front, which might have different settings than the cliff critical profile. It is not possible to obtain the depth across the profile from the images or, for example, estimate the length of the cracks on the back of the cliff overhang. Figure 10 shows the maximum tensile stresses as function of a range of soil depths (0-2m). The estimates are calculated for Profiles 3-1 to 3-5. The cliff section tensile stresses increase significantly due to the increase of the soil cover thickness on the top of the cliff (Fig. 10, profiles 3-3 and 3-4). For example, an increase of 160% in the tensile stresses is obtained if the soil depth increases from 0 (no soil) to 1m. Hence, further work should focus on improvements to these estimates. There are some limitations to this methodology, which will require further consideration in the future. An assumption was made that the density is constant along the whole cliff face, despite breccias being lithologically heterogeneous rocks and therefore having a spatially varied density. The density can also vary with time. Table 3. The effect of the notch length on estimates of the tensile stresses. Profile 5 Property Symbol Unit ln_min ln_max ln_mean Notch depth ln m 5.02 6.56 5.79 Cliff height hc m 6.92 6.92 6.92 Bulk density Y kN/m3 25.30 25.30 25.30 Momentum M kNm/m1 2206 3767 2935 Centre of cliff masse Z m3 7.98 7.98 7.98 Tensile stress (equations 1-3) o kN/m2/m1 276.40 472.00 367.70 Tensile stress (real area) o kN/m2/m1 284.30 284.30 284.30 Tensile stress difference % -3% 66% 29% So far the stability method, adopted here, does not take into account other parameters and effects that might influence the stability of cliffs. Secondary processes such as rock weathering, soil erosion and subaerial processes are not included. Furthermore, the method does not take into account processes that might influence the formation and stability of notches, such as wave impact, salt-water weathering, wetting and drying and beach sediment abrasion. Figure 7 illustrates that the sea level reaches the notches at several locations. The narrow gravel beach absorbs most of the wave's mechanical energy during normal conditions. However, during an extremely high water level, waves reach the cliff's toe and can contribute to the cliff's undercut formation and the cliff's retreat. According to the expected sea-level rise and more frequent events of extreme high water levels (Aqua Alta) in the Kvarner area (Antonioli et al., 2007), the sea level and waves will reach the cliff toe more frequently. This could contribute to the acceleration of marine erosion and cliff retreat. It is important to investigate how to incorporate the dynamic processes of changes in cliff geometry into the stability model. 6 CONCLUSIONS The investigated area around the Stara Baška settlement is in a delicate geodynamic balance. The comparison of ortho-rectified aerial photos from 1966 and 2004 showed that marine erosion is quite prominent and the cliff retreat has been between 4 and 5 metres with a simultaneous expansion of the beach surface. The main cause of the cliff slopes' instability is the occurrence of strong waves and the formation of wave-cut notches. Secondary causes are the weathering and erosion of the soil and rocks on the slope surface of the cliff. The slump of the cliff's slope can occur in a rock mass with better strength parameters, where the notches are cut a few metres inward into the toe of the cliff slope. The new combined method for the stability analysis of coastal cliffs was tested and had incorporated the cantilever-beam model and the structure-from-motion (SfM) photogrammetry. The SfM photogrammetry can provide highly detailed 3D cliff geometrical data, which can then be utilized with the stability model. This is particularly important in a stability analysis of rocks such as breccias with a varied geometry, which cannot be easily replaced by a rectangular surface. The SfM method overcomes the limitations of traditional surveys in estimating the cliff overhang surface and the notch length. It has been demonstrated that the tensile stresses estimated by using 'real' and 'assumed rectangular' cliff areas may differ significantly. Acknowledgements This paper is a result of a research activity under the project "Geological hazard in the Kvarner area" by the University of Rijeka. REFERENCES [1] Abam, T.K.S. 1997. Genesis of channel bank overhangs in the Niger Delta and analysis of mechanisms of failure. Geomorphology 18, 2, 151-164. [2] Antonioli, F., Silenzi, S. 2007. Variazioni relative del livello del mare e vulnerability delle pianure costiere italiane. Quaderni della Societa Geologica Italiana 2, 1-29 (in Italian). [3] Benac, Č., Juračić, M. 1998. Geomorphological indicators of sea level changes during upper Pleistocene (Würm) and Holocene in the Kvarner region (NE Adriatic Sea). Acta Geographica Croatica 33, 27-45. [4] Benac, Č., Juračić, M., Bakran-Petricioli, T. 2004. Submerged tidal notches in the Rijeka Bay NE Adriatic Sea: Indicators of relative sea-level change and of recent tectonic movements. Marine Geology 212, 21-33. [5] Benac, Č., Juračić, M., Matičec, D., Ružić, I., Pikelj, K. 2013. Fluviokarst and classical karst: Examples from the Dinarics (Krk Island, Northern Adriatic, Profile 3 1 2 3 4 5 Subprofile Number Figure 10. Different soil layer (on the top of the cliff, Fig. 3) height and cliff critical profile tensile stresses (profiles 3-1 to 3-5). Croatia). Geomorphology 184, 64-73. [6] Booij, N., Ris, R.C., Holthuijsen, L.H. 1999. A third-generation wave model for coastal regions, Part I: Model description and validation. Journal of Geophysical Research 104 (C4), 7649-7666. [7] Castedo, R., Murphy, W., Lawrence, J., Parede, C. 2012. A new process-response coastal recession model of soft rock cliffs. Geomorphology 178, 128-143. [8] Davies, P., Williams, A.T., Bomboe, P. 1991. Numerical modelling of Lower Lias rock failures in the coastal cliffs of South Wales. Coastal Sediments '91. - American Society of Civil Engineers, New York, pp. 1599-1612. [9] Griggs, G.B., Trenhaile, A.S. 1997. Coastal cliffs and platforms. Coastal Evolution: Late Quaternary shoreline morphodynamics, Cambridge, pp. 425-450. [10] Hoek, E., Brown, E.T. 1991. Practical Estimates of Rock Mass Strength. International Journal of Rock Mechanics and Mining Science 34, 8, 1165-1187. [11] Hoek, E. 1996. Strength of Rock and Rock Masses. ISRM News Journal 2, 2, 4-16. [12] Hungr, O., Evans, S.G., Bovis, M., Hutchinson, J.N. 2001. Review of the classification of landslides of the flow type. Environmental and Engineering Geoscience 7, 221-238. [13] James, M.R., Robson, S. 2012. Straightforward reconstruction of 3D surfaces and topography with a camera: Accuracy and geoscience application. Journal of Geophysical Research - Earth Surface 117(F3), F03017. doi:10.1029/2011JF002289. [14] James, M.R., Ilić, S., Ružić, I. 2013. Measuring 3D coastal change with a digital camera. Proc. of the Coastal Dynamics 2013. - 7th International Conference on Coastal Dynamics, Arcachon, pp. 893-904. [15] Juračić, M., Benac, Č., Pikelj, K., Ilić, S. 2009. Comparison of the vulnerability of limestone (karst) and siliciclastic coasts (example from the Kvarner area, NE Adriatic, Croatia). Geomorphology 107, 1-2, 90-99. [16] Kogure, T., Aoki, H., Maekade, A., Hirose, T., Matsukura, Y. 2006. Effect of the development of notches and tension cracks on instability of limestone coastal cliffs in the Ryukyus, Japan. Geomor-phology 80, 3-4, 236-244. [17] Kogure, T., Matsukura, Y. 2010. Critical notch depths for failure of coastal limestone cliffs: case study at Kuro-shima Island, Okinawa, Japan. Earth Surface Processes and Landforms 35, 9, 1044-1056. [18] Kuhn, D., Prüfer, S. 2014. Coastal cliff monitoring and analysis of mass wasting processes with the application of terrestrial laser scanning: A case study of Rügen, Germany. Geomorphology 213, 153-165. [19] Leder, N., Smirčić, A., Vilibić, I. 1998. Extreme values of surface wave heights in the Northern Adriatic. Geofizika 15, 1, 1-13. [20] Marinos, P., Hoek, E. 2000. GSI: A geologically friendly tool for rock mass strength estimation. Proceedings of the GeoEng 2000, Melbourne, pp. 1422-1446. [21] Matsukura, Y. 1988. Cliff instability in pumice flow deposits due to notch formation on the Asama mountain slope, Japan. Zeischrift für Geomorpolo-gie 32, 129-141. [22] Matsukura, Y. 2001. Rockfall at Toyohama Tunnel, Japan, in 1996: effect of notch growth on instability of a coastal cliff. Bulletin of Engineering Geology and the Environment 60, 285-289. [23] Pikelj, K., Juračić, M. 2013. Eastern Adriatic Coast (EAC): Geomorphology and Coastal Vulnerability of a Karstic Coast. Journal of Coastal Research 29, 4, 944-957. [24] Podobnikar, T. 2009. Methods for visual quality assessment of a digital terrain model. S.A.P.I.EN.S [Online], 2, 2, URL: http://sapiens.revues.org/738 [accessed on April 4th 2014] [25] Ružić, I. 2003. Analysis of the Sea level interaction with the rivers high water levels in the North Adriatic, the Dubračina mouth example. BSc thesis, University of Rijeka, Faculty of Civil Engineering, Rijeka, (in Croatian). [26] Ružić, I., Marović, I., Vivoda, M., Dugonjić Jonjčević, S., Kalajžić, D., Benac, Č., Ožanić, N. 2013. Application of Structure-from-Motion photogrammetry for erosion processes monitoring, Moscenicka Draga example. Proc. of the 4th Workshop of the Japanese-Croatian Project on "Risk Identification and Land-Use Planning for Disaster Mitigation of Landslides and Floods in Croatia, University of Split, Split, pp. 49-50. [27] Ružić. I., Marović, I., Benac, Č., Ilić, S. (2014). Coastal cliff geometry derived from structure-from-motion photogrammetry at Stara Baška, Krk Island, Croatia. Geo-Marine Letters. doi 10.1007/s00367-014-0380-4. [28] Snavely, N., Seitz, S.N., Szeliski, R. 2008. Modeling the world from internet photo collections. International Journal of Computer Vision 80, 189-210. [29] Sunamura, T. (1992). Geomorphology of Rocky Coasts. Wiley, Chichester. [30] Thorne, C.R., Tovey, N.K. 1981. Stability of composite river banks. Earth Surface Processes and Land-forms 6, 469-484. [31] Timoshenko, S.P., Gere, J.M. 1972. Mechanics of Materials. Van Nostrand Reinhold Co., New York. [32] Trenhaile, A.S. 2002. Rock coasts, with particular emphasis on shore platforms. Geomorphology 48, 1-3, 7-22. [33] Westoby, M.J., Brasington, J., Glasser, N.F., Hambrey, M.J., Reynolds, J.M. 2012. "Structure-from-Motion" photogrammetry: A low-cost, effective tool for geosci-ence applications. Geomorphology 179, 300-314. [34] Williams, A.T., Davies, P., Bomboe, P. 1993. Geometrical simulation studies of coastal cliff failures in Liassic strata, South Wales, UK. Earth Surface Processes and Landforms 18, 703-720. NOVA METODA ZA PREIZKUŠANJE TRDNOSTI GLINE NA PORUŠITEV ZARADI PRECE-JANJA PRI VISOKIH VODNIH TLAKIH Izvleček Porušitev gline pri visokih vodnih tlakih je težko oceniti. У članku je predstavljena nova metoda za oceno porušitve gline s trdnostjo proti precejanju, ki predstavlja kritični vodni tlak, da razruši glino. Načrtovan je bil eksperiment za hitrejše preizkušanje vrednosti precejne deformacije gline glede na tradicionalno metodo. У članku sta predstavljena sistem in metoda eksperimenta, prikazani so rezultati 18 skupin vzorcev, ki so bili preizkušeni. Rezultati kažejo, da na trdnost proti precejanju vplivata debelina glinenega sloja in precejne poti. Pokazalo se je tudi, da se voda infiltrira v glino v pogojih, ko vodni tlak preseže minimalno vrednost (P0). Fu-wei Jiang (vodilni avtor) Guizhou Institute of Technology, School of Resources and Environmental Engineering Kitajska E-pošta: jfwei_666@126.com Ming-tang Lei Institute of Karst Geology and Key Laboratory of KrastCollapse Prevention, CAGS, Kitajska E-pošta: mingtanglei@hotmail.com Xiao-zhen Jiang Institute of Karst Geology and Key Laboratory of KrastCollapse Prevention, CAGS, Kitajska E-pošta: jiangxiaozhen2005@hotmail.com Ključne besede porušitev gline; precejna deformacija; trdnost (upor) proti precejanju; visok vodni tlak A NEW METHOD FOR TESTING THE ANTI-PERMEABILITY STRENGTH OF CLAY FAILURE UNDER A HIGH WATER PRESSURE Keywords clay failure; seepage deformation; anti-permeability strength; high water pressure Abstract It is difficult to judge the failure of clay seepage under a high water pressure.This paper presents a new method to assess clay failure based on the anti-permeability strength, which is the critical water pressure to destroy the clay. An experiment is designed to test the value that avoids the problem of the time-consuming, traditional method to test clay seepage deformation. The experimental system and the process of testing are introduced in this paper. With a self-designed experimental system and method, 18 groups of sample were tested. The results show that the clay thickness and the seepage paths influence the anti-permeability strength. It also indicates that water infiltrates into the clay under the condition that its pressure exceeds a minimum value (P0). 1 INTRODUCTION Clay is an important structural material to reduce hydraulic conductivity and can be used widely for liners and caps in many projects. It is usually safe to use as an anti-seepage material, but some cases of clay failure have occurred. At present, the research on clay failure concentrates on the hydraulic conductivity (K) and the surface erosion. The stability of clay projects is usually judged by the value of K [1-6]. K is a parameter for which the appropriate values are difficult to obtain with good accuracy due to a large pore system [7-8], the effect of the measuring scale [9], the condition including freezing-thawing and drying-wetting [10], overlying soil [1], chemical composition [11], water content [12], etc. The K values measured with an oedometer are lower than those in a triaxial test and vary a great deal with the applied pressure [10], and reduce with an increasing density of the sample [6]. Surface erosion is mainly a form of clay failure. In nature, clay erosion accompanies the evolution of landforms. Couper [13] demonstrated how the surface processes vary with the silt-clay content of river-bank soil and how to consider this variation in the context of erosion observed in the field. Lévy [14] discussed how erosion and landslides are related to valley development. Lamelas [15] shares some of the characteristics of previous models, including the erosion of bare clay surfaces by wave-generated, bottom shear stresses and of mobile, sediment-covered surfaces by abrasion. Another model demonstrates that bluff height, debris mobility, wave undercutting, and groundwater levels are key factors in determining the clay stability of coastal bluffs [16]. In addition, clay erosion influences cultivation. Clay particles can be eroded from the plough layer and Fu-wei Jiang (corresponding author) Guizhou Institute of Technology, School of Resources and Environmental Engineering China E-mail: jfwei_666@126.com Ming-tang Lei Institute of Karst Geology and Key Laboratory of KrastCollapse Prevention, CAGS, China E-mail: mingtanglei@hotmail.com Xiao-zhen Jiang Institute of Karst Geology and Key Laboratory of KrastCollapse Prevention, CAGS, China E-mail: jiangxiaozhen2005@hotmail.com transported both laterally and vertically, through pores and cracks into the backfill, and then directly to drain pipes [17]. Messing [9] carried out an in-situ experiment to quantify the effect of a plough on the displacement of soil down aslope. Its result shows that the clay soil removed from the top of the plot is much greater than that estimated for the surface-water erosion. Mukonen [18] advises improving the soil structure and avoiding the disruption of stabilized aggregates in order to prevent erosion in an agriculture field. However, clay failure includes some forms, except surface erosion, one of which is internal erosion under a high water pressure. Though a dam elevates the level of the water and produces high pressures, the dam core is often made of clay so as to prevent any seepage. However, this has potential risks [19-22]. Previous research on the issue of a clay core is mainly concerned with the material character [23-25], stability [26], testing pressure [27], contact boundary [28], stress strength [29], erosion and permeability [30-31], and monitoring [32]. The process of internal clay erosion, however, is known to contribute to dam failure during high water pressure.The high pressure can lead to an excessive ingress of water into the clay layer, ending in breaching and failure. But very little information is known about internal erosion and how to quantitatively evaluate it. There are no established methods to solve the problem. The seepage-failure criterion is usually the critical hydraulic gradient. But it is difficult to judge the clay seepage failure, because it has been proved by laboratory soil tests that the seepage deformation of clay is nonlinear. So the critical hydraulic gradient remains unfitted. The time for that water to penetrate into the clay is the seepage-failure criterion of the clay in this study. Generally, the coefficient of permeability of the clay is less than 10-6cm/s. Based on the equation of Darcy's law of V=KI, it shows V is less than 10-3cm/s when 7=1000, and it means that water takes at least 1000sec to penetrate a 1-cm thickness of clay under the condition of K=10-6cm/s. It represents a clay failure that the water takes a long time to penetrate normally into the clay. That is to say, the clay structure is a failure when water takes seconds or minutes to penetrate into the clay. The penetrating time is related to water pressure, such that the higher the water pressure, the quicker the water penetrates into the clay. By increasing the water pressure, when water penetrates through the clay at the bottom in seconds or in minutes, the pressure is critical, which is the anti-permeability strength. This study regards the time as the critical condition to judge the clay's failure and introduces a quick method to test its value quantitatively. 2 EXPERIMENTAL APPARATUS AND METHOD 2.1 Experimental apparatus The experimental apparatus (Fig.1) is designed to test the anti-permeability strength of clay under a high water pressure. The system (Fig.2) consists of a triaxial instrument, a water pipe, valves, and a sample pot. The confining pressure of the triaxial instrument is chosen to support a high water pressure. The instrument is a Model TCK-1 Triaxial Test Controlled Apparatus, which is made by the Nanjing Soil Instrument factory, one of the earliest and largest enterprises specializing in the production of geotechnical engineering laboratory instruments in China. It performs very well in a high-pressure cell composed of a confining chamber (confining pressure up to 2.5 MPa), and is convenient for reading the value via a digital display. The water conduit connects the triaxial instrument to the sample pot. The water pressure is transmitted along the water conduit from the confining chamber to the soil Model TCK- 1 Triaxial Test Controlled Apparatus Water Conduit Air Outlet Valve Small Stone Pricked Hole Inlet Valve Drainage Outlet Sample Pot Soil Sample Figure 1. Schematic diagram of the experimental apparatus. triaxial instrument Д л Figure 2. Testing equipment. sample. In order to regulate effectively the water pressure, it installs the inlet valve to control the transmission of pressure by opening or closing the valve. The sample pot consists of a steel cover, a plastic pipe (diameter 91mm) and rubber packing. They work via a screw. At the top of the plastic pipe, an air outlet valve is fixed to release the air in the water conduit and the sample pot when injecting water into the test equipment. The center tail of the sample pot produces a hole, with a diameter of 10mm, to observe whether the water seeps and the particle is out. It is only under sealed conditions that the triaxial instrument can support the confining pressure of the water. In the design system, the clay sample at the end of the testing system is most likely opening, but it is not for two reasons. Clay can stop the water loss for the hardly impermeable character. Air is removed in the whole process of the testing program. There is no water and air loss, so it is impossible to release the water pressure unless the clay fails. 2.2 Experimental Procedures In this study, the procedure for the experiment has seven steps, as follows: a) Sample preparation. The sample is placed in the sample pot, and the contact surfaces are smeared with Vaseline to prevent water or air leakage. b) Hole piercing.The conical-headed piercing device (20 cm long; 1.5 mm diameter) (Fig. 3) is used to pierce ♦ Figure 3. Pricking device. _*_ Figure 4. Pricked hole. a hole in the center of the sample (Fig. 4). This allows the water to flow, just as water seeps naturally via weaknesses in a subterranean environment. c) Installation. Following the approach shown in Fig. 1 and the described method, the installation is tested to ensure functionality. d) Saturation. The sample is saturated using a standard vacuum-pump method. e) Air purging. Firstly, open the inlet and exhaust valve of the sample pot. Secondly, fill with water until the air is completely purged from the test equipment. Finally, close the valve. f) Increasing water pressure. Rotate the knob controlling the confining pressure at a constant velocity to increase the water pressure (1KPa/s). The water pressure is up to 1 MPa in a few minutes. It is important to be careful as this step is the key to accomplishing the test successfully. h) Observation and recording. With increasing pressure, observe the phenomenon of water seepage, such as the permeating capacity of the water and the water's turbidity, bubbles, and the flowing particle, and record the corresponding pressure value, crucially. 3 JUDGING THE CONDITION OF ANTI-PERMEABILITY STRENGTH The water pressure provided by the confining chamber of the Model TCK-1 Triaxial Test Controlled Apparatus increases quickly (1 KPa/s), even if the knob is rotated slowly. It works on clay and enlarges the pressure difference between the pores that results in the pores expand- ing and the structure of the clay particle changes until failure. It has the minimum water pressure that results in the clay's failure, which is judged by the clay particle flows or the coefficient of permeability becomes larger. The former is difficult to observe. The latter appears when water seeps through a certain thickness of clay in a short time as the clay has a low coefficient of permeability (<10-9m/s) and an impermeable layer. Fig.5 shows the condition for judging the clay's failure by a variation of the permeability coefficient. In addition, the water pressure is in a closed environment. Once it is destroyed, the water pressure decreases to zero. As the sample is a part of the testing system, that it fails results in the water pressure reducing at the original constant speed of the increasing pressure. The pressure at the time when the water pressure starts to decrease in testing process is the anti-permeability strength. 4 MATERIAL TESTED The specimens tested in this study come from the industrial park of Guilin, China, N 25°13'17.51", before Figure 5. The condition of judging clay failure. E 110°14'31.50". Its natural physical property given in Table 1 is tested by the soil laboratory of the Institute of Karst Geology. The clay deposited during the Quaternary is characterized by a hard plastic shape, a block structure, and a high porosity.The average number of the standard penetration test is 9 on this layer, with a thickness of 7m. The water level lies at adepth of 3-4m. As the clay lying water level is weak compared to the other layer, the testing program utilizes the clay layer at a depth of 3.5 meters in order to achieve its purpose. Actually, the sample of undisturbed clay is unsuitable for testing due to the high porosity that results in a deceasing water pressure. It compresses the clay to improve the tested sample. The water content (w), the dry density (pd), and the void ratio (e) of the compacted sample are 31.27%, 1.635g/m3 and 0.671, respectively. 5 RESULTS Table 2 shows the testing results for the18 experiments with horizontal and vertical types, 1.0, 1.5, 2.0 and 2.5cm Table 2. Testing result of anti-permeability strength. Sample No. Seepage type Thickness/cm Anti-permeability strength/KPa Average value/KPa Y1 62 Y2 1.0 68 64.3 Y3 63 Y4 73 Y5 1.5 66 68.0 Y6 n о N ir 65 Y7 О Я 79 Y8 2.0 85 79.7 Y9 75 Y10 72 Y11 2.5 86 80.0 Y12 82 Y13 85 Y14 ^ 1.5 92 85.0 Y15 c 78 Y16 ■e 87 Y17 2.0 90 90.7 Y18 95 Table 1. Physical property of original clay. w % p kg/m3 pd kg/m3 Gs kg/m3 e n % Sr % WL % WP % Ip Il C MPa ф° Es MPa 36.1 1810 1330 2740 1.064 51.6 93.0 48.0 28.6 19.4 0.38 0.050 19.3 7.22 of thickness. All the tested samples were pricked with a 0.5-cm-deep hole (the red shown in Fig.1) at the center with a pricking deviceto lead the infiltrating way, for the hole is the weakest area that is easily permeated.The testing results show that the seepage types and the sample thickness influence the anti-permeability strength. The anti-permeability strength is increasing with thickness in the horizontal test type (Fig.6) as well as in the vertical test type (Fig.7). The relationship liner in Fig.6 is matched well with the average value for R2=0.8891. It Original Value Average Value -linear (Average Value) i I- 65 60 55 Po JO У = 11,73 R- = Зх— 52.4 3.8891 67 / ► / * ( \ / к t У' К < У Л / с i 0,0 0,5 1,5 2,0 The thickness of sample (cm) Figure 6. Anti-permeability strength of different thicknesses with the horizontal test type. does not start from the original point, for the function has aconstant (52.467) as the initial pressure (P0) caused by some factors. One factor is the bonding water film around the clay grains. The water film prevents the water from permeating into the internal channel. Thus, the initial pressure reflects the resistance of the water film. When the water pressure over comes the resistance, the water permeates into the clay. The other factor is the particle force, including the friction, the molecular gravity, and the electrostatic force. Clay has a strong particle force to prevent the water from permeating strongly, which is why clay has a lower permeability. Water should exceed the particle force firstly, and then permeate into the clay. Therefore, P0 is the minimum pressure to make the water infiltrate into the clay. The anti-permeability strength varies with different thicknesses (Fig.7). In nature, the seepage deformation of the clay varies. The horizontal and vertical types that are common are stimulated in this paper. The vertical type only chooses the downward route. Compared to the strength of the horizontal type (Table 2), the vertical average strengths of the 1.5 and 2.0 cm samples are increased by 25%, and 13.8%, respectively. They are not the same because of the different P0 caused by the different resistances. The resistance comes from the gravity. In the downward direction, the gravity promotes water seepage. But, the gravity resists water seepage in the horizontal type. It can be predicted that the strength of the downward direction is less than the horizontal type, as its resistance to gravity is small. The critical hydraulic gradient (i=H/L, Fig.8) is nonlinear based on the results data. It indicates that the water seeping in clay media exhibits non-Darcian flow behavior, which is the same result as other researchers [33-36]. So the critical hydraulic gradient does not fit to judge * Horizontal test type ■ Vatical test type 100 ss 95 t 90 - 1 85 !> i SO i 75 Cu 1 70 65 60 1 1,5 2 1 The thickness ofsample(cm) Figure 7. Anti-permeability strength of different thickness. Figure 8. Critical hydraulic gradient. the clay failure. There are four factors that explain the problem, as follows: a) the pores of the internal clay are very small, which results in a strong capillary force that influences the effect of the water head; b) in the process of clay-seepage deformation, the water pressure drives the water to infiltrate into the clay while it compresses the sample that decreases the porosity and changes the structure of the sample; c) the influence of the bonding water film and the particle force was mentioned above; d) the distribution of the effect stress caused by the water pressure is heterogeneous, and it is unsuitable to apply the theory of Darcy's Law in this case. The water pressure makes the clay fail at its weakness point, not the whole, and the hydraulic gradient is the volume force. Therefore, it is difficult to obtain the effect value that the hydraulic gradient should not be applied to all the samples. 6 DISCUSSION Water seepage in clay media is governed by Richards' equation [37], with the proportionality factor termed the hydraulic conductivity (K). It is both difficult and time consuming to measure the functional relationship between K and the soil water-pressure head (H) [7, 8, 38]. Experimental results indicate that the critical hydraulic gradient is not suitable directly to judge clay failure. The main reason is that the traditional form of Darcy's law is not appropriate for describing water-flow processes in clay media, because the observed relationship between the water flux and the hydraulic gradient can be very non-linear [39]. The non-Darcian flow behaviors in low-permeability media were also investigated by Klausner [40], Alabi [41, 42], while the complexity of field-scale water flow in a clay formation was discussed in [43]. Clay failure depends not only on the hydraulic gradient, but on the friction, electrostatic force, molecular attraction and capillary water pressure between the pores. The critical hydraulic gradient described the relation between the seepage force and gravity, and does not reflect the complex distribution of stress caused by the water pressure [44]. The physical indicators of clay, such as the permeability coefficient, the critical hydraulic gradient, etc., changes while the clay is compressed and the degree of consolidation increases under a high water pressure [45]. Therefore, it is inappropriate to simply consider that the clay failure is determined by the hydraulic gradient. In this paper, the anti-permeability strength was applied to judge the clay failure under a high water pressure. Compared to the hydraulic gradient, it comprehensively reflects the mechanical properties of clay seepage for it responds to the influence of structure, compression, particle force, etc. It makes up for the disadvantage of the traditional measure, to some extent. Moreover, the values of the experimental results are between 20 and 80 KPa, which is equal to the 2-8m of water head that is common in Nature. Consequently, it is effective and positive to judge clay failure by its anti-permeability strength. Although the method is not perfect, it gives a new view and a quick way to assess the risk of clay failure under a high water pressure. The experimental apparatus is easy to install in most geotechnical laboratories. Experimentalists can master the testing process quickly. The water pressure is up to hundreds of kilopascal (KPa) in minutes using this system. That makes up for the deficiency of the traditional method to test clay-seepage deformation. It is significant to apply the method of testing anti-aging strength to assess the stability of a dam's foundation made with clay, the probability of reservoir leakage, the susceptibility of clay cave collapse in a Karst region, the pollution risk in a reservoir area due to leakage, etc. However, the mechanism of anti-permeability strength needs to be researched and the correlation between the strength and the sample size, including the ratio of the height to the diameter, the size of the cross-sectional area, the porosity, the saturation level, etc., should be explored further. 7 CONCLUSION The experimental apparatus and a quick method are designed carefully to test the anti-permeability strength of clay under a high water pressure in this study. Model TCK-1 Triaxial Test Controlled Apparatus providing a high water pressure in minutes in the system avoids the time-consuming problem, which is the weakness of the traditional method to test clay-seepage deformation. The method gives an effective new view to assess the risk of clay failure under a high water pressure. A total of 18 groups of samples were tested with the self-design test equipment and test method. It interprets the relationship between the anti-permeability strength and the clay thickness and seepage ways. The results show that: a) the anti-permeability strength is increasing with the thickness; b), it exists at the minimum pressure (P0) to make water infiltrate into the clay; c), compared the horizontal type, the vertical average strength of the 1.5 and 2.0cm samples are increased by 25%and 13.8%, respectively. Acknowledgement This work is supported by the National Natural Science Foundation of China (No.41302288). We also sincerely thank You-qiang QIN for providing the natural physical property of the tested material. REFERENCES [1] Fox, P. J., De Battista, D.J., Mast, D. G. 2000. Hydraulic performance of geosynthetic clay liners under gravel cover soils. Geotextiles and Geomem-branes 18, 2, 179-201. [2] Bhardwaj, A.K., Goldstein, D., Azenkot, A., Levy, G.J. 2007. Irrigation with treated wastewater under two different irrigation methods: Effects on hydraulic conductivity of a clay soil. Geoderma 140, 1, 199-206. [3] Sui, W.H., Liu, J. Y., Du,Y. 2009. Permeability and seepage stability of coal-reject and clay mix. Proce-dia Earth and Planetary Science1, 1, 888-894. [4] Shackelford, C.D., Sevick, G. W., Eykholt, G.R. 2010. Hydraulic conductivity of geosynthetic clay liners to tailings impoundment solutions. Geotex-tiles and Geomembranes 28, 2,149-162. [5] Cho, W.J., Lee, J.O., Choi, H. J. 2012. Radionuclide migration through an unsaturated clay buffer under thermal and hydraulic gradients for a nuclear waste repository. Annals of Nuclear Energy 50, 71-81. [6] Hamdi, N., Srasra, E. 2013. Hydraulic conductivity study of compacted clay soils used as landfill liners for an acidic waste. Waste Management 33, 1,6066. [7] Bouma, J., Wösten, J.H.M. 1979. Flow patterns during extended saturated flow in two, undisturbed swelling clay soils with different macro-structures. Soil Sci. Sot. Am. J. 43, 16-22. [8] Germann, P., Beven, K. 1981.Water flow in soil macropores. I. An experimental approach. J. Soil Sci. 32, l-13. [9] Messing, I., Jarvis, N.J. 1995. A comparison of near-saturated hydraulic properties measured in small cores and large monoliths in a clay soil. Soil Technology 7, 4, 291-302. [10] Graham, J., Yuen, K., Goh, T. B., Janzen, P., Siva-kumar, V. 2001. Hydraulic conductivity and pore fluid chemistry in artificially weathered plastic clay. Engineering Geology 60,1, 69-81. [11] Zhang, X.C., Norton, L.D. 2002. Effect of exchangeable Mg on saturated hydraulic conductivity, disag-gregation and clay dispersion of disturbed soils. Journal of Hydrology 260, 4, 194-205. [12] Berilgen, S.A., Berilgen, M.M., Ozaydin, I.K. 2006. Compression and permeability relationships in high water content clays. Applied Clay Science 31, 3, 249-261. [13] Couper, P. 2003. Effects of silt-clay content on the susceptibility of river banks to subaerial erosion. Geomorphology 56, 2, 95-108. [14] Lévy, S., Jaboyedoff, M.,Locat, J., Demers, D. 2012. Erosion and channel change as factors of landslides and valley formation in Champlain Sea Clays: The Chacoura River, Quebec, Canada. Geomorphology 145,12-18. [15] Lamelas, M. T., Hoppe, A., de la Riva, J., Marinoni, O. 2009. Modelling environmental variables for geohazards and georesources assessment to support sustainable land-use decisions in Zaragoza (Spain). Geomorphology 111, 1, 88-103. [16] Castedo, R., Fernandez, M., Trenhaile, A.S., Pare-des, C. 2013. Modeling cyclic recession of cohesive clay coasts: Effects of wave erosion and bluff stability. Marine Geology 335, 162-176. [17] 0ygarden, L., Kvsrner, J., Jenssen, P. D. 1997. Soil erosion via preferential flow to drainage systems in clay soils. Geoderma 76, 1, 65-86. [18] Mukonen, P., Hartikainen, H., Alakukku, L. 2009. Effect of soil structure disturbance on erosion and phosphorus losses from Finnish clay soil. Soil and Tillage Research 103, 1, 84-91. [19] Vaughan, P.R., Soares, H.F. 1982. Design of filters for clay cores of dams. International Journal of Rock Mechanics and Mining Sciences &Geome-chanics Abstracts 19, 4, 88. [20] Capozio, N.U.O., Dupuls, M.M. 1983. Geotechni-cal problems related to the building of a tailings dam on sensitive varved clay. International Journal of Rock Mechanics and Mining Sciences &Geome-chanics Abstracts 20, 5, 156. [21] Nagarkar, P.K.K., Kulkarni, M.V., Kulkarui, D.G. 1984. Failures of a monozone earth dam of expansive clay. International Journal of Rock Mechanics and Mining Sciences &Geomechanics Abstracts 21, 3, 103. [22] Wang, W.S. 1986. Earthquake damages to earth dams and levees in relation to soil liquefaction and weakness in soft clays. International Journal of Rock Mechanics and Mining Sciences &Geome-chanics Abstracts 21, 3, 118. [23] Brookins, D.G. 1972. Possible accumulation of authigenic, expandable-type clay minerals in the substructure of tuttle creek dam, Kansas, U.S.A. Engineering Geology 6, 4, 251-259. [24] Dingsov, G., Markov, G., Alexieva, L. 1975. Estimation of the consolidation pattern of waterproof clay cores of earth dams from local materials. Inter- national Journal of Rock Mechanics and Mining Sciences &Geomechanics Abstracts 12, 4, 62. [25] Rengasamy, P., McLeod, A.J., Ragusa, S.R. 1996. Effects of dispersible soil clay and algae on seepage prevention from small dams. Agricultural Water Management 29, 2, 117-127. [26] Fry, J. J, Delage, P., Nedjat, N., Nanda 1993. Computing the stability of clay fill dams under construction. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 30, 5, 313. [27] Charles, J. A., Watts, K. S. 1987. Measurement and significance of horizontal earth pressures in the puddle clay cores of old earth dams. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 24, 3, 115. [28] Stematiu, D., Popescu, R., Luca, E. 1991. Contact clay problems during the erection of Maru dam. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 28, 5, 318. [29] Valore, C. 1991. Strength parameters of a tectonized clay for the design of a large dam. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 28, 2, 138. [30] Atkinson, J.H., Charles, J.A., Mhach, H.K. 1990. Examination of erosion resistance of clays in embankment dams. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 27, 6, 338. [31] Muijs, J.A., Kruse, G.A.M. 1991. Erosion and permeability of material for clay liners on dikes. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 28, 6, 352. [32] ul Haq, I. 1993. Monitoring of clay core: Mangla Dam project. International Journal of Rock Mechanics and Mining Sciences&Geomechanics Abstracts 30, 5, 313. [33] Hansbo, S. 1960. Consolidation of clay, with special reference to influence of vertical sand drains. In: Swed. Geotech. Inst. Proc., 8. [34] Swartzendruber, D. 1961. Modification of Darcy's law for the flow of water in soils. Soil Sci. 93, 22-29. [35] Miller, R.J., Low, P.F. 1963. Threshold gradient for water flow in clay systems. Soil. Sci. Soc. Am. Proc. 27, 6, 605-609. [36] Blecker, R.F. 1970. Saturated Flow of Water Through Clay Loam Subsoil Material of the Brol-liat and Springerville Soil Series. The University of Arizona, MasterThesis. [37] Richards, L.A. 1931. Capillary conduction of liquids in porous mediums. Physics, 1, 318-333. [38] Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Res. 12, 513-522. [39] Liu, H.H., Birkholzer, J. 2012. On the relationship between water flux and hydraulic gradient for unsaturated and saturated clay. Journal of Hydrology 475, 242-247. [40] Klausner, Y., Craft, R. 1966. A capillary model for non-Darcy flow through porous media. J. Rheol. 10, 603-613. [41] Popoola, O.I., Adegoke, J.A., Alabi, O.O., 2009. Modification of fluid flow equation in saturated porous media. Global J. Pure and Appl. Sci. 15, 395-400. [42] Alabi, O.O. 2011. Validity of Darcy's law in laminar regime. Electron. J. Geotech. Eng. 16, 27-40. [43] Antia, D.D.J. 2008. Oil polymerization and fluid expulsion from low temperature, low maturity, overpressured sediments. J. Pet. Geol, 31, 3, 263-282. [44] Zhou, H.X., Cao H. 2011. Research on mechanism and numerical simulation method of seepage of double-layer dike foundation. Chinese Journal of Rock Mechanics and Engineering 10, 2128-2136. [45] Zhuang, X.S., Zhao, X., Zhu, R.G. 2005. Researching cohesive soil critical hydraulic gradient change law under the action of Penetration.Urban Geotechnical Investigation & Surveying 3, 45-47. NUMERICNO MODELIRANJE NEGATIVNEGA TRENJA PLAŠČA NA POSAMIČNEM VERTIKALNEM IN NAGNJENEM PILOTU Izvleček У članku je preučevano negativno trenje na posamičnem vertikalnem in nagnjenem pilotu. Razvit je bil model s končnimi elementi (z uporabo ABAQUS programske opreme) na osnovi študij Lee in drugi in Comodromos in Bareka. Rezultati modela so bili primerjani z znanimi potrjenimi rešitvami, nato je bil za različne obtežbe na površini analiziran posamični vertikalni pilot, ki nosi na nogi in plašču, pri čemer sta bili preučevani trenjska sila zaradi negativnega trenja in položaj nevtralne ravnine. Nadalje je bila izvedena analiza odziva posamičnega pilota, ki nosi na nogi in plašču za različne kote naklona pilota med 0 in 30°. Poleg tega je bila z 2-D modeli izvedena analiza občutljivosti. Analiza je pokazala zadovoljivo skladnost z rezultati študije avtorjev Lee in drugi ter Comodromos in Bareka. Na koncu je bilo ugotovljeno, da sta trenjska sila zaradi negativnega trenja in položaj nevtralne ravnine odvisna od stanja tal okoli pilota, od uporabe 2-D ali 3-D modela, intenzitete obtežbe na površini, tipa pilota (nosilnost na nogi ali na plašču) in odklona pilota od vertikale. Mohammad Mahdi Hajitaheri Ha Imam Khomeini International University, Civil Engineering Department Qazvin, Iran E-pošta: mahdi_hajitaheri@yahoo.com Mahmoud Hassanlourad Imam Khomeini International University, Civil Engineering Department Qazvin, Iran E-pošta: hassanlou@eng.ikiu.ac.ir Ključne besede negativno trenje na plašču, nagnjen pilot, trenjska sila zaradi negativnega trenja, nevtralna ravnina NUMERICAL MODELING OF THE NEGATIVE SKIN FRICTION ON SINGLE VERTICAL AND BATTER PILE Abstract In this paper the negative skin friction on single vertical and batter pile is investigated. First, a finite-element model (using ABAQUS software) based on the studies Lee et al. and Comodromos and Bareka was developed. After that, the results of the model were compared and validated. Then a single vertical end bearing and a single skin friction pile under different earth-surface loadings were analyzed and the down-drag force as well as the neutral plane location were studied. Subsequently, the performances of the single end bearing and the friction pile, with different inclination angles between 0 and 30°, were analyzed. Moreover, the sensitivity analysis was implemented using 2-D models. This showed a satisfactory compatibility with the results of the study of Lee et al. and Comodromos and Bareka. Finally, it was concluded that the drag load of the pile and the neutral plane position depend on the condition of the soil surrounding the pile, the 2D or 3D model type, the earth-surface loading intensity, the type of end-bearing pile or friction pile and the pile's inclination angle. The simulation results agree well with the experimental findings. Mohammad Mahdi Hajitaheri Ha Imam Khomeini International University, Civil Engineering Department Qazvin, Iran E-mail: mahdi_hajitaheri@yahoo.com Mahmoud Hassanlourad Imam Khomeini International University, Civil Engineering Department Qazvin, Iran E-mail: hassanlou@eng.ikiu.ac.ir Keywords negative skin friction, batter pile, down-drag force, neutral plane 1 INTRODUCTION Evaluation of Negative Skin Friction (NSF) is a common problem in the design and construction of a pile in soft ground that is subjected to an earth-surface loading or the drawdown of an underground water table, etc. [1, 2, 3]. The development of an additional compressive force (dragload) in the form of excessive pile settlement (downdrag) could cause construction and maintenance difficulties for the supported structure [4]. However, the NSF is a function of the normal effective stress on the pile skin and for the development of friction forces, two surfaces are needed to prevent a contact between the two surfaces, but the drag load cannot cause settlement by itself [5]. The NSF is usually mobilized at its ultimate limit, except for the proximity to the neutral plane, and a limiting relative displacement of the soil and the pile is assumed to achieve full mobilization of the skin friction [6, 7]. So, disregarding the design recommendations of the pile drag load could lead to a collapse of the structure, e.g., a quay or bridge, built on pile foundations. Several reports can be found relating to the NSF, the dragload and the down drag on single piles. Poulos and Davis (1980) [8], Chow et al. (1990) [9], Lee (1993) [10], Teh and Wong (1995) [11], Fellenius (1972 and 1989) [12, 13] indicated that neutral planes that are developed in piles depend on the surrounding soil's settlement and the rigidity of the pile material; they also mentioned that the larger the bearing of the pile tip, the deeper the neutral plane, and accordingly, the higher the drag force. Jinyuan et al. (2012) [14] studied the NSF in a single pile under various influencing factors, including the consolidation time, the pile/soil interface, the lateral earth-pressure coefficient, the soil-pile limiting displacement, the intensity of the surcharge , the soil compressibility, and the stiffness of the bearing layer. It was found that the neutral plane changes significantly with the consolidation time and the stiffness change in the bearing layer. Van Der Veen (1986) [15] stated for a given pile, the more allowable the settlement, the deeper the neutral plane, and accordingly, the lower the drag force. The majority of these reports have introduced semi-empirical formulae to predict the location of the neutral plane, and accordingly to estimate the drag force for a given pile and soil condition. However, these formulae do not incorporate the settlement of the pile as a governing parameter. However, the drag load of batter piles has not been studied so widely. The batter pile is applied in combination or exclusively with a vertical pile as the transmission of force to the foundation to be used. The batter pile's performance is a suitable action against the horizontal loads. Very limited information and experimental data can be found about the batter piles' NSF in the literature. Hanna and Nguyen (2003) [17] showed that the shaft resistance decreases with an increase in the inclination angle of the pile for both compression and tension loadings. The full-scale test results reported by Meyerhof and Yalsin (1992) [16] showed that the shaft resistance capacity of a pile increases due to an increase in the pile's inclination angle, while Hanna and Afram (1986) [26] concluded that there is no significant change in the shaft's resistance with an increase in the pile's inclination angle. The results showed that the effect of the type of pile and the load distribution were an important parameter that could affect the bearing capacity and the pile's resistance. Lee et al. (2002) [18] studied the drag-load, the down-drag and discussed the efficiency of pile groups under NSF conditions using a FEM, and the effects of parameters like the ground surface load, the friction coefficient of the pile-soil interface, the arrangement of the pile groups, and the spacing of the piles. Emilios et al. (2005) [19] studied the drag-load and the neutral plane of a single pile in a layered soil influenced by the sequence of the pile-head and the earth-surface loading using FLAC3D software. Poulos (2008) [20] presented a design method for pile foundation under NSF, and discussed the NSF's effective characteristics on the residual stress of the pile, the pile-head loading and the efficiency of the pile group. Prakash and Subraman-yam (1965) [21] developed an experimental batter pile model and indicated that the resistance of the negative batter pile is more than that of the positive batter pile in the lateral loading condition. Rajashree and Sitharam (2001) [22] described the nonlinear modeling of vertical and batter piles according to the finite-element method. In this model, the resistance and the lateral displacement of the batter pile induced by static and cyclic loads were investigated. In this research the soil was not modeled as continuous media, but was considered as springs with equivalent stiffness. Poulos and Madhav (1971) [23] considered the manner of a batter pile in a lateral static load and showed that laterally the displacement of a batter pile depends on the pile's inclination angle. In this research some experimental results of Emilios et al. (2005) [19] were compared with the finite-element model obtained using ABAQUS software. First, a 3D numerical model was used to determine the negative skin friction and the drag load in a vertical friction and end-bearing pile under different earth-surface loadings. Then some effective parameters, like the loading and the inclination angle of the pile, were studied and compared in 2D and 3D conditions. The soft clay was undrained in the numerical model. 2 MODELING AND CALIBRATION The effects of negative skin friction on the drag load and the down drag were studied and the batter and vertical piles were compered in the modeled soil according to Lee et al. (2002) [18]. The pile was considered to be elastic with a diameter of 0.5 m and a height of 20 m. The profile of the soil was layered as a top soft clay layer that was underlined by a bottom dense sand layer. The pile shaft was placed in the soft soil and the end of it was settled on the dense sand. The elements that were used in the three-dimensional model of the soil and the pile are continuous eight-node brick elements, the reduced integration method (C3D8R) and continuous 20-node cubic elements (C3D20). In addition, continuous 4-node plain stress elements (CPS4) were considered in the two-dimensional model. The constitutive models that were used for the soil and the pile are the elastic-perfectly plastic Mohr-Coulomb non-associated flow rule model and the linear elastic model, respectively. A fine mesh was employed for the central part of the soil, where the maximum stresses and strains are expected. Following a short parametric study, a total of 8 brick elements around the cylinder circumference in this central part were found to be sufficient to achieve a convergence of solution, whereas the size of the brick elements in the longitudinal direction was chosen to be equal to 0. 5 m. The soil and pile materials' properties are summarized in Table 1. Geometry, stratification and meshes of the model made using the 3D finite-element software ABAQUS is shown in Fig 1. Table 1. Material properties used in the analysis by Lee et al. (2002). Material Model E (kPa) C (kPa) V ф (0) w (0) K0 Y (kN/m3) Concrete pile Isotropic elastic 2,000,000 0.15 25 Soft clay Mohr-Coulomb 5000 30 0.4 20 0 0.65 18 Dense sand Mohr-Coulomb 50,000 0.1 0.3 45 10 0.5 20 20 m 5 m Elastic behaviour Ycrit Displacement Figure 1. Finite-element model of a vertical single pile in soil. Figure 2. Behavior of the interface element by Lee et al. (2002). The size of the soil elements in the vicinity of the pile is finer to increase the accuracy of the analysis, and then larger further from the pile to avoid a long computational time. The boundary of the model is far enough from the pile to resist the stress conflict. The model was first solved to achieve the in-situ stresses and then the initial displacements were set to zero and before the main analysis, the pile and soil interaction was used with the penalty method and the interface friction coefficient was supposed to be 0.3-0.4. Also, the soil's residual stresses change due to pile's installation was overlooked. Because of the type of pile, according to the method of installation, it was assumed to be bored. The stiffness of the lower soil was assumed to be the same as the upper soft-soil layer to model the frictional pile in the NSF analysis. However, it is assumed to be 10 times larger than the upper soft layer for the end-bearing pile. The constitutive model of the interface elements between the pile and soil is defined by the Mohr-Coulomb shear failure criterion. The interface elements can be separated if tension is developed across the interface and gapping is formed between the pile and the soil interface and the shear and normal forces are set to zero. The normal effective stress, a', between two contact surfaces was multiplied by the interface friction coefficient, p, to give the limiting frictional shear stresses т = p •a'. If the shear stress applied along the interface is less than т = p •a, the surfaces will be stuck .Otherwise, the nodes of the elements in the contact will slide (Fig.2). It was found that the pile's behavior is governed mainly by the interface behavior. The coefficient of the Earth Pressure at Rest (K0) was considered using the method of Jaky (1948) [25] for normal soils for clay and dense sand as 0.65 and 0.5, respectively. The interface friction coefficient, p, (p = tan S) was supposed to be 0.3-0.4, because it was accepted in the Mohr-Coulomb way. As shown in Figs. 3 and 4, the shear stress on the pile skin and the axial drag load of the pile resulting from the present numerical model are very close to Comodromos (2005) [19] and Lee. (2001) [24]. Simplified methods assumed that whole mobilization of the shear strength above and below the pile neutral line; however, an analysis using ABAQUS software let us consider the partial mobilization of the shear strength according to the applied shear strain. S-L = 50 kpa Dragload (KN) 50 kPa -200 -400 -600 Friction pile, Lee. et al. (2002). Present Model \ / a) Shear Stress (kPa) 150 100 50 0 -50 ' -100 -150 I I I I I I I I I I I I I I I fciO I I I I I I I I I I I I I d e D - 20 -Frictionpile, Lee); et al. (2002). - Present Model ■ Neutral Plane L 25 b) Figure 3. a) Shear-stress distributions along the pile surface and b) axial drag load for a single friction pile. S-L = 50 kPa 0 200 Dragload (KN) 400 600 800 -End-Bearing pile. Lee. et al. (2002). — End-Bearing Pile, Comodromos, (2005). ---3D (present model) S-L = 50 kPa 0 10 0 5 ) e 20 - 25 -1 Shear Stress (kPa) 20 30 40 50 I i i i i I i i i i I Friction pile, Lee. et al. (2002) a) b) Figure 4. a) Shear-stress distributions along the pile surface and b) axial drag load for a single end-bearing pile. 3 COMPARISON OF 2D AND 3D MODELS In general, the combination of the axial load and the moment would cause the pile to behave differently than the vertically loaded pile; therefore, the critical case for this combination was dependent on the type of piles. In forming the plastic inside batter pile we tend to bend the plastic to be created within the pile behavior batter pile proportionally different. 2D and 3D modeling of the vertical and the end-bearing batter piles were first compared under an earth-surface loading of Ap = 50kpa. As is clear from Fig. 5, there is compatibility between the results of the 2D and 3D modeling and the overall trend of skin friction along the pile is similar. The skin shear stress at the pile head is equal to zero in the 2D modeling, which differs from the 3D modeling. The 2D modeling shows less skin shear stress up to a pile depth of about 5 m; however, it shows more skin shear stress from a depth of 5 m up to nearly the pile end, in comparison with the 3D modeling. Both types of modeling show approximately the same results at the pile end. The influence of strain in the third dimension has mainly the effect in the case of soil failure. This factor S-L = 50 kPa Shear Stress (kPa) 0 10 20 30 40 50 60 O I I I I I I I I I I I I I I I I I I I I I I I I I I 5 - I 10 h t a e О 1? 20 25 -L Figure 5. Skin shear-stress development along the end-bearing vertical pile. Dragload (KN) 2000 25 —Ж—Skpa —*— 50 KPa ---150 KPa 25 KPa Figure 6. Drag load along the single friction pile under various earth-surface loads. could affect the actual bearing pile capacity. In the 3D modeling, at the end of the pile shear, the failure zone spreads out in three directions; neverthless, the shear zone is two dimensional in the 2D modeling. This means a larger end bearing capacity and less skin friction development for the pile in the 3D modeling, compared to the 2D modeling. 4 DRAG LOAD AND NEUTRAL PLANE IN THE VERTICAL SINGLE PILE In this part, the magnitude of the drag load and the location of the neutral plain are estimated for the friction and the end bearing pile under different earth-surface loads. As shown in Fig. 6, increasing the earth-surface load leads to an increase in the cumulative drag load. The application of a surface load of 5 KPa leads to a maximum drag load of 300KN, and increasing the surface load up to 200 KPa shows an increasing trend in the drag load up to a maximum value of 1600KN, indicating a non-linear relation between the surface load and resulting drag load. The decreasing trend after the maximum value of the drag load is due to the change of the shear-stress direction (the neutral line) in the pile skin. Since in this case the neural point in the end-bearing pile is formed lower and subsequently the cumulative force in the shell of pile against the down drag is less, it leads to more bearing capacity. Fig.7 presents the displacement of the neutral plane along a friction pile for an earth-surface loading of 5 kPa S-L = 5 kPa 0 Settlement (mm) 2 4 -........; a <г л / Ì : Ü —©— Soil settlement Pile settlement a) SL = 200 kPa 0 0 5 ) (m h £ 15 D 20 25 Settlement (mm) 50 100 150 i i i i i i ........t 1 / ® * я / J 9 a ^ J / CT I __a__cpttlpmpnt Pile settlement b) Figure 7. Pile and surrounding soil settlements for a surface load of: a) 5 and b) 200 KPa. and 200 kPa. It is clear that the location of the neutral plane changes from a depth of 14m to 17.4m, as the surcharge increases from 5 KPa to 200 KPa. Fig. 8 shows the drag load imposed on the end-bearing vertical pile for various earth-surface loadings. 5 DRAG LOAD IN THE BATTER PILE In order to investigate the effect of the pile's inclination drag load a numerical model was produced in 3D conditions. The inclination angle of modeled batter pile was 15 degrees with respect to the vertical axis. The 3D model is shown in Fig 9. Figure 8. Drag load along the single end bearing pile for different surface loads. Figure 9. Finite-element mesh of a single batter pile and a soil profile. As expected, larger drag loads are developed as earth-surface surcharge increases. As the pile is rested on a stiff layer with a surrounding soft soil, the neutral plane forms near the stiff layer so that an increased surcharge can push the neural plane down to the end of the pile. Comparing Figures 6 and 8 indicated that under a constant earth surcharge, the drag load is more for an end-bearing pile than for a friction pile. The material properties, including the soil and the pile, completely correspond to those mentioned in Fig.1 and Table 1. It is clear that the boundaries were defined far enough from the pile in order to eliminate the interaction between the boundaries and the pile. The pile was assumed to be loaded at the rest condition (k0), before inducing the surcharge. As the surcharge is (b) Mohr-coulomb Circule c^cl o' Figure 10. a) A cross-section of a soil profile batter pile and the applied stress on the pile, b) limits of the normal pressure acting on he pile skin. loaded, the batter pile is subjected between the active (ka) and the passive (kp) loading, which produces a complex condition to predict the behavior of the pile. It is possibly related to the flexural deflections and the longitudinal distortions of the pile. Fig. 10 illustrates the explained process for a batter pile. Based on this figure, between ka and kp along the pile for a batter pile. 6 EFFECT OF THE INCLINATION ANGLE IN THE DRAG LOAD FORCE_ For a closer look at the batter pile behavior under different inclination angles, some sensitivity analyses for both the friction and the end-bearing batter piles under a constant earth-surface load of Ap = 50kPa were performed. The loading of the batter pile is similar on the vertical pile loading for comparing the results under the same conditions. Figure 11 shows the drag force of the end-bearing piles for different inclination angles. The sign and direction of the drag force change along the pile and this condition intensifies with an increase of the pile inclination angle. In this condition the pile is subjected between the active and the passive force surrounding from the soil, the skin shear stress and then the drag force's direction change along the pile. In other words, the pile is subjected to a longitudinal distortion and flexural bending that may lead to a change in the sign of the shear stress and the drag force in it. The deformed shape of the batter pile is shown in Fig.12, illustrating the above process. This is the reason why the batter pile is subjected to a different circumferential force at the different depths. However, it is certain that a greater inclination angle leads to a greater dragload of the batter pile. Figure 12. Deformation of an end-bearing batter pile under earth-surface loading. Figure 13 shows the drag force of the batter frictional pile with different inclination angles. As shown, the pile drag force increases with an increase of the inclination angle. It is considered that the maximum drag load force occurs nearly at the end of the frictional batter pile and then rapidly decreases. The friction batter pile's behavior is simpler with respect to the sign of the drag force and the pile's longitudinal distortion. Comparing Figures 11 and 13 we can conclude that the drag force of the frictional batter pile is about four times larger than that of the end-bearing batter pile. S-L = 50 kPa -200 0 -ft Dragload (KN) 200 400 600 25 Figure 11. Drag load along the end-bearing batter pile. S-L = 50 kPa 0 5 Dragload (MN) 10 15 20 25 Pile angle 5° —©—10» —й—20° ---30* 25 L Figure 13. Drag load along the friction batter pile. 7 CONCLUSIONS REFERENCES This paper deals with the negative skin friction in a single vertical and a batter pile in two-layer soil. A model based on the study of Lee et al. (2002) [18] was developed and the results were compared and evaluated. The results of the model showed good compatibility with this study's conclusions. Then some other factors, such as the third-dimension of the model, the intensity of the surface load, the batter pile's inclination angle, and the type of pile were investigated. In this study the soil conditions were in accordance with the paper of Lee et al. (2002), so the soil characteristics can vary as a factor in future discussions. The analysis of the results shows that: 1. The drag load of the pile is completely dependent on the surrounding soil's condition, the intensity of the surface load and the batter pile's inclination angle. Also, the neutral plane's position on the pile depends on these conditions and its location tilts into the resistant subsoil layer. 2. In the vertical or near-vertical positions, the pile acts as a column in the soil and is only subjected to an axial force However, for large inclination angles, close to 25 degrees, the pile works like a beam column that is subjected to an axial force, the flexural momentum sand shear forces. In this situation the normal forces between the pile and the soil vary between the active and passive cases and make the interaction between the pile and the soil more complex. 3. The drag load of the pile under an earth-surface loading for an end-bearing vertical pile is more than a pile. This trend is also dominant in the batter pile; the difference is that the maximum drag load in the pile does not necessarily occur in the neutral plane. 4. Increasing the earth-surface loading, the neutral plane in the frictional and the end-bearing pile moves to a greater depth of the soil. The neutral plane is almost at the bottom of the pile for the end-bearing piles; however, it is formed in the end one-third height of the piles for the frictional pile. 5. Despite the vertical pile, no specific neutral plane forms in the batter piles. In the batter piles, the shear-stress direction on the pile periphery changes in a region. Also, the pile drag load increases with an increase of the inclination of pile. Acknowledgements The authors want to acknowledge Dr. M. Mousavi, assistant professor of K.N. Toosi University of Technology, who edited this paper [1] Phamvan, P. 1989. Negative skin friction on driven piles in Bangkok subsoils. PhD thesis, AIT, Bangkok, Thailand. [2] Little, J.A., Ibrahim, K.I. 1993. Predictions associated with the pile downdrag study at the SERC soft clay site at Bothkennar, Scotland. In Predictive Soil Mechanics (Wroth Memorial Sympo- sium), London: Thomas Telford 796-818. [3] Lee, C.J., Chen, H.T., Wang, W.H. 1998. Negative skin friction on a pile due to excessive groundwa-ter with drawal Proc. International Conference, Centrifuge 98, New York, 513-518. [4] Fellenius, B.H. 1999. Piling terminology. http:// www.geoforum. com/info/pile info/terminology. asp. [5] Fellenius B.H. 1984. Negative Skin Friction and Settlement of Piles. Second International Seminar, Pile foundations, Nanyang Technological Institue, Singapore. [6] Broms, B. 1979. Negative skin friction. Proc. 6th Asian Regional Conf. Soil Mech. Found. Engng, Singapore 2, 4-75. [7] Fellenius, B.H. 2006. Results from long-term measurement in piles of drag load and downdrag. Can. Geotech Journal 43(4), 409-430. [8] Poulos, H.G., Davis, E. H. 1980. Pile foundation analysis and design, Wiley, New York. [9] Chow, Y.K., Chin, J. T., Lee, S. L. 1990. Negative skin friction on pile groupsInt. J. Numer. Analyt. Meth. Geomech., 14, 75-91. [10] Lee, C.Y. 1993. Pile groups under negative skin friction. J. Geotech. Eng. 119(10), 1587-1600. [11] Teh, C.I., Wong, K.S. 1995. Analysis of downdrag on pile groups Geotechnique, 45(2), 19-207. [12] Fellenius, B.H. 1972. Down-drag on piles in clays due to negative skin friction. Can. Geotech. J. 9(4), 325-337. [13] Fellenius, B.H. 1989. Unified design of piles and pile groups Transportation Research Record 1169, Transportation Research Board, Washington, D.C, 75-82. [14] Jinyuan, L., Hongmei, G., Hanlong, L. 2012. Finite element analyses of negative skin friction on a single pile, Springer J. Acta Geotechnica 7, 239-252. [15] Van Der Veen 1986. A general formula to determine the allowable pile bearing capacity in case of negative friction DFI, Proc., Int. Conf. on Deep Foundations, Beijing, China, 2,138-147. [16] Meyerhof, G.G., Yalsin, A.S 1992. Behavior of flexible batter Pile under inclined loads in layered Soil. Can. Geotechnical J. 30, 247-256. [17] Hanna, A., Nguyen, T.G. 2003. Shaft Resistance of Single Vertical and Batter Piles Driven in Sand. ASCE J Geotechnical and Geoenviromental Eng.129, 601-607. [18] Lee, C.J., Bolton, M.D., Al-Tabbaa, A. 2002. Numerical modeling of group effects on the distribution of drag loadsin pile foundations, J. Geotechnique 52 (5), 325-335. [19] Comodromos, E.M., Bareka, S. V. 2005. Evaluation of negative skin friction effects in pile foundations using 3D nonlinear analysis. Computers and Geotechnics 32, 210-221. [20] Poulos, H.G. 2008. A partial design approach for pile with negative friction. GEI Geotechnical Eng. Issue 161(2), 19-27. [21] Prakash, S., Subramanyam, G. 1965. Behavior of battered piles under lateral loads. J. Indian Nat. Soc of Soil Mech. And Found. Eng., New Delhi, 4,177196. [22] Rajashree, S. S., Sitharam, T.G. 2001. Nonlinear finite element method of batter piles under lateral load. J. Geotech. Eng. 127(7), 604-612. [23] Poulos, H.G., Madhav, M.R. 1971. Analysis of movement of battered piles. Rep. No. R173, University of Sydney, 1-18. [24] Lee, C.J. 2001. The influence of negative skin friction on piles and in pile groups. Ph.D. thesis, Cambridge University. [25] Jaky, J. 1948. Pressure in silos. Proc., 2nd Int. Conf. on Soil Mechanics and Foundation Engineering, Rotterdam, the Netherland, 1, 103-107. [26] Hanna, A.M., Afram. A. 1986. Pull-out capacity of single batter piles in sand. Canadian Geotechnical Journal 23(3), 387-392. NASIP NA MEHKIH TLEH OJAČAN S PILOTI ZA ŽELEZNICE VISOKIH HITROSTI: NUMERIČNA IN ANALITIČNA RAZISKAVA Izvlecek Z geosintetiki ojačan pilotni nasip predstavlja učinkovito in gospodarno metodo za reševanje težav glede nosilnosti temeljnih tal in nesprejemljivih posedkov ter nestabilnosti brežin nasipov zgrajenih na mehkih tleh, kar vodi do široke uporabe, zlasti za nasipe železnic visokih hitrosti. Razvite so bile nekatere metode projektiranja za ocenjevanje uspešnosti teh ojačanih konstrukcij, ki v glavnem temeljijo na rezultatih iz modelov v zmanjšanem merilu in numeričnih simulacij. Zanesljivost teh metod mora biti preverjena s terenskimi preizkusi v naravni velikosti. V članku je predstavljena numerična in analitična študija na terenskem preizkusu v naravni velikosti, in sicer, na nasipu za železnico visokih hitrosti Fengyang. Rezultati so bili analizirani in komentirani za posedke zemljine, razmerje koncentracije napetosti (RKN), osno silo in trenjsko napetost v pilotu. Pokazalo se je, da so rezultati posedkov dobljeni z metodo končnih elementov (MKE) in analitično metodo bili zelo blizu merjenim, in to je bil zanesljiv parameter, da bo ocenitev obnašanja s piloti ojačanega nasipa sprejemljivo natančna. RKN je bil po modificirani Terzaghijevi metodi precenjen za 25%, medtem, ko je bil po MKE podcenjen za približno 20%. Pokazano je bilo tudi, da bi se natezna sila v ojačitvi lahko učinkovito ocenila s predlagano analitično metodo, medtem ko je bila s FEM precenjena za 44%. Yan Zhuang (vodilni avtor) Hohai University, Geotechnical Research Institute Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering 1 Xi Kang Road, Nanjing210098, Kitajska E-pošta: zhuangyan4444@hotmail.com Xiaoyan Cui Hohai University, College of Civil and Transportation Engineering 1 Xi Kang Road, Nanjing210098, Kitajska E-pošta: cui19890213@126.com Ključne besede nasip ojačan s piloti, železnica visokih hitrosti, numerična simulacija, analitična metoda REINFORCED PILED EMBANKMENT FOR A HIGH-SPEED RAILWAY OVER SOFT SOIL: A NUMERICAL AND ANALYTICAL INVESTIGATION Keywords reinforced piled embankment, high-speed railway, numerical simulation, analytical method Abstract A geosynthetic, reinforced, piled embankment is an effective and economic method to solve the problems of possible bearing failure, unacceptable settlement and slope instability for an embankment built over soft soil; this has led to its widespread use, especially for high-speed railway embankments. Some design methods have been developed to assess the performance of these reinforced structures, which are mainly based on the results from small-scale models and numerical simulations. However, the reliability of these methods needs to be validated under full-scale field tests. This paper presents a numerical and analytical study for a full-scale field test of the Fengyang high-speed railway embankment. The results were analyzed and discussed in terms of the settlement of subsoil, the stress-concentration ratio (SCR), the axial force and the fric-tional stress of the pile. They showed that the settlement of the subsoil, from both the finite-element method (FEM) and the analytical method, were in good agreement with the measurement, and thus was a reliable parameter to assess the performance of the piled embankment with reasonable accuracy. The SCR was overestimated by the modified Terzaghi method, with a difference of 25%, while it was underestimated by the FEM, with a difference of approximately 20%. It was also shown that the tensile force in the reinforcement could be effectively assessed using the proposed analytical method, while it was overestimated by the FEM with a difference of 44%. Yan Zhuang (corresponding author) Hohai University, Geotechnical Research Institute Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering 1 Xi Kang Road, Nanjing210098, China E-mail: zhuangyan4444@hotmail.com Xiaoyan Cui Hohai University, College of Civil and Transportation Engineering 1 Xi Kang Road, Nanjing210098, China E-mail: cui19890213@126.com 1 INTRODUCTION The high-speed railway is the result of modern science and technology, and an important indicator of railway modernization. The construction of a high-speed railway network in China has been developing rapidly over the past 10 years, for example, the Beijing-Shanghai, Shanghai-Hangzhou, Shanghai-Nanjing, and Wuhan-Guangzhou high-speed railways. The controls of the post-construction settlement and the improvement in construction speed have been the two main technical problems for high-speed railways constructed on soft soil. The piled embankment has been proven to be an interesting solution that prevents the failure or excessive deformation of the embankment compared to the traditional soft-foundation improvement methods for embankments built over soft soils, such as preloading (Magnan, 1994)[1]. The inclusion of a geosynthetic reinforcement enhances the load-transfer mechanism and so minimizes the differential and maximum settlements (Gangakhedkar, 2004) [2]. Many experimental tests, criteria and analytical methods have been proposed to analyze the interactions between the embankment, the geosynthetic reinforcement, the piles and the subsoil. Terzaghi (1943)[3] assumed that the shearing resistance of the soil during arching was mobilized along two vertical planes through the side of the inclusion, and the considered shear stress on the vertical interfaces originating from the rigid supports on either side of a 'trapdoor, where the support was reduced. Russell and Pierpoint (1997) [4] developed the arching theory proposed by Terzaghi (1943) [3] to take account of the three-dimensional problems of the piled embankments. Hewlett and Randolph (1988) [5] conducted 3D model tests and presented a semi-circular model (in 2D) or a hemispherical dome model (in 3D) to describe the arching, and of the uniform thickness, with no overlap. This method assumed that the pressure acting on the subsoil was uniform and considered the failure condition either at the crown of the arch or the pile cap. The British Standard BS8006-1 (2010) [6] estimated the vertical stress on the top of the piles based on the modified Marston's formula for positively projecting conduits. The magnitudes of the load carried by the piles and the reinforcement were evaluated. The three-dimensional shape analysis developed by Hewlett and Randolph (1988) [5] was also included in BS8006-1 (2010) [6]. However, the effect of the subsoil in the two methods of BS8006-1 (2010) [6] was not considered. The criterion for the German EBGEO (2011) [7] was based on the method proposed by Kempfert et al. (2004) [8]. In this method, the average vertical pressure acting on the columns and the soft foundation soil were calculated by considering hemispherically shaped domes spanning the distance between the pile caps. Van Eekelen et al. (2011) [9] analyzed and modified BS8006-1 (2010) [6] in terms of the calculation method of the tensile force in the tensile reinforcement for the 3D situation. They reported that BS8006-1 (2010) [6] designed a relatively strong and expensive geosynthetic reinforcement in comparison with other design models (e.g., EBGEO, (2011) [7]), and the modified method was better in evaluating the geosynthetic reinforcement than that of the original BS8006-1 (2010) [6], but the geosynthetic reinforcement predicted by EBGEO (2011) [7] was much closer to the FEM results. Zhuang et al. (2014) [10] proposed a simple analytical model for the reinforced-piled embankment to assess the contribution of the reinforcement and the subsoil. The simplified model can be used to estimate the magnitude of the tensile force generated in the reinforcement, the vertical stress increment acting on the subsoil for multi-layered soft soils, the maximum settlement at the surface of the subsoil and the maximum strain of the geosynthetics. The methods illustrated above are either based on the results from small-scale models or the numerical analysis, whose reliability should be validated under full-scale field tests. Some research has also been conducted based on a numerical simulation, which has the advantage of being able to analyze the interactions between the embankment, the geosynthetic reinforcement, the piles and the soft-soil foundation with reasonable accuracy. Han and Gabr (2002) [11] performed a numerical study of reinforced, piled embankments, including the underlying subsoil. Le Hello and Villard (2009) [12] introduced a numerical and experimental method for a geosynthetic reinforced, piled embankment. Two main mechanisms, including the arching effect and the membrane effect of the geosynthetics, were analyzed. But very limited attention has been paid to a settlement-time relationship, especially a post-construction settlement, which is critical to the performance of pavements, especially a high-speed railway built on embankments. This paper presents a numerical and analytical study for the field case of the Fengyang high-speed railway. The results are discussed in terms of the stress-concentration ratio (SCR), the development of the settlement of subsoil at the centerline of the embankment, the distribution of the settlement of subsoil along the width of the embankment, the lateral deformation of the subsoil, the axial force and the negative frictional stress along the piles. The purpose of this paper is to compare the results of the numerical simulation and the analytical methods with measurements made on a full-scale field case to assess their accuracy in analyzing the behavior of the reinforced, piled embankment. 2 DESCRIPTIONS OF THE FENGYANG HIGHSPEED RAILWAY CASE STUDY 2.1 General project information The test site of the Fengyang high-speed railway was located in Fengyang County, Anhui Province, China. Details of its construction, the ground conditions, the setup of the tests, the instruments and the measurement results can be obtained from Zeng (2010) [13]. In the present paper, only the information required to describe the numerical simulation and the analytical methods was proposed. The embankment is 3.6-4.8 m high with a crest width of 13.6 m. The side slope of the embankment is 1V to 1.5H. The thickness of the soft soil ranges from 13.1 to 16.6 m. Below the soft soil is the weathered rock, with a bearing capacity in the range 250 to 800 kPa. The CFG (cement fly-ash gravel) pile-net composite foundation was used in the field case to improve the soft ground. The geosyn-thetic reinforcement with a tensile stiffness of 1 MN/m in the longitudinal and transverse directions was also used to strengthen the embankment. The field case of the Fengyang high-speed railway was tested with a total length of 370 m from DK854+640 to DK855+010. The embankment was divided into three sections, with different ground-improvement methods. In this paper, only the section of DK854+655 (defined as Section A), reinforced by CFG pile-net composite foundation and preloading, was analyzed. The preloading in Section A was about 3 m high with a side slope of 1V to 1H, and was applied to the top of the embankment for approximately 450 days. Detailed information about Section A is summarized in Table 1. 2.2 Site conditions The height of the embankment in Section A was 3.6 m. The fill material consisted mainly of gravel mixed with sand, with an average unit weight of 20 kN/m3. The CFG piles were arranged in a square pattern with a spacing of 1.8 m, and designed with circle pile caps (see Fig. 1(a)). Table 1. General information of Section A for the Fengyang high-speed railway. Section A Embankment height (m) Pile length (m) Pile spacing (m) Pile cap diameter (m) Cushion structure 3.6 10.0 1.8 1.0 0.6 m gravel cushion + a layer of geogrid ©Inclinometer П Settlement template A Earth pressure cell J) Strain sensor .__13-6 m c) 1.0 Settlement template \ 5 Cushion b) Strain sensor Figure 1. Plan view and cross-section of the test embankment in the Fengyang high-speed railway: a) Plan view of the experimental site at Section A; b) Cross-section of the embankment; c) Structure of the cushion. Table 1 summarizes the information about the piles. Fig. 1(b) shows the cross-section of the test embankment. It can be seen that two layers of clay were laid below the embankment. The bearing capacity of the first layer of clay was 150 kPa, with a thickness of 6.5 m, while the thickness of the second layer was 9.0 m, with a bearing capacity of 200 kPa. The main properties of the clays are summarized in Table 2. Beneath the clay layers are two layers of weathered rock with a high bearing capacity. One layer of the biaxial geogrid was sandwiched between two, 0.3-m-thick, gravel layers to form a 0.6-m-thick composite-reinforced bearing layer (see Fig. 1(c)). Table 2. Main soil properties of Section A. Stra- Water con- Void Compres- Critical shear tum tent (%) ratio sion index stress ratio Clay 1 25.0 0.8 0.23 0.79 Clay 2 22.7 0.7 0.19 1.03 2.3 Test setup and instruments As shown in Fig. 1, the earth pressure cells were placed on the top of the pile caps and the surface of the subsoil. Piezometers to measure the pore-water pressures were placed at different depths of the subsoil. The settlement plates were installed on the top of both the pile caps and on the surrounding soils. Strain sensors in the geogrid were installed to measure the strain of the geogrid. The observation started at the beginning of the construction for the embankment in November in 2007 and lasted for about 2 years after it was completed. The measurements of the field case will be presented and discussed in the following sections to validate the results derived from the numerical simulation and analytical methods, in terms of the vertical stress acting on the pile cap and the subsoil, the settlement and the lateral deformation of the soil, the axial force and the negative frictional stress along the piles. 3 NUMERICAL SIMULATION_ 3.1 Description of the numerical model The finite-element method (FEM), which was carried out using the finite-element package 'ABAQUS' (version 6.12), was used to simulate the field case of the Fengyang high-speed railway. Since the embankment was symmetrical along its centerline, only half of the embankment was modeled for the FEM. The foundation soil was taken to be 21.5 m deep, overlying a rigid impermeable stratum. The horizontal length of the model was taken to be 36.6 m, so that the boundary effect can be minimized. A 1.8-m-wide section with a center row of piles beneath the embankment was selected for the analysis (see Fig. 2(a)), so that a truly full-3D model can be obtained. The 3D numerical model is shown in Fig. 2(b). The water table was at a depth of 0.8 m below ground level, and the initial pore pressures prior to the embankment construction were taken to be hydrostatic. A zero-pore-pressure boundary condition was applied at the level of the water table to model the free drainage, which means that the water is allowed to drain via this plane of the soft clay. The bottom of the model is defined as impermeable, and a lateral flow is not permitted across the vertical planes. With regard to the displacement boundary condition, no displacements in the direction perpendicular to the symmetry planes and to the base were allowed. The embankment fill was assumed to be dry, granular material, and was modeled with a linear-elastic, perfectly 3 О о о о о о I/' роооооо ) о о о о о о u a) 1 SmV Figure 2. Three-dimensional finite-element model of the embankment: (a) Plan view for layout of CFG piles; (b) 3-dimen-sional finite-element model. Table 3. Summary of the material parameters used in the finite-element analyses. Material E (MPa) Y (kN/m3) c (kPa) ф (0) v Embankment 60 20 40.0 20.0 0.15 Clay 1 6.4 19 37.3 20.4 0.30 Clay 2 9.2 20 75.6 27.3 0.30 Weathered rock 1 40 21 63.9 19.3 0.25 Weathered rock 2 80 22 63.9 19.3 0.25 CFG pile 22000 23 0.10 Geosyntheic reinforcement Tensile stiffness kg = 1.0 MN/m; Poisson's ratio v = 0.3 (assumed) plastic model with the Mohr-Coulomb failure criterion, as well as two layers of weathered rocks. The pile and the geogrid were modeled as a linear elastic material. The two clay layers were modeled using a modified cam clay (MCC) model. The parameters used are summarized in Tables 2 and 3. The critical shear stress ratio of the two clay layers (see Table 2) for the MCC model were 0.79 and 1.03, respectively, and with a Poisson's ratio value of 0.3. The compressibility of the subsoil can be quantified by the compression index (Cc). In cam clay, the compressibility is represented by the parameters X and к. The parameter X is equal to Cc divided by 2.3, and it is assumed in this study that к is equal to 0.1 times X. In addition, the initial yield surface size a0, a required input for the MCC model, is computed using the following: 1 ao = 2exp (1- eN - eoklnpo I — (1) where e0 is the initial void ratio, p0 is the overburden pressure, eN is the intercept of the normally consolidated line with a ratio axis in the e-lnp' plane, X is the slope of the virgin consolidation line, and к is the slope of the swelling line. The geogrid in the FEM was modeled using 4-node, full-integration, 3-dimensional membrane elements (M3D4). The clay layers below the water table were modeled using 8-node, stress-pore-pressure-coupled, brick elements (C3D8P), while 8-node brick elements (C3D8) were used for the other materials. The interface friction angle between the geogrid and the embankment was assumed to be equal to the friction angle of the embankment fill, which was also assumed in Liu et al. (2007) [14]. According to the equation proposed by Potyondy (1961) [15] for the clay, the interface friction angle (ф) between the pile and the subsoil can be determined as follows: j = (0.6~0.8) j (2) where fs is the friction angle of the soil. In this analysis, the ratio of the interface friction angle (ф) and the friction angle of the soil (ф5) was assumed to be 0.7. And duplicated nodes were used to form a zero-thickness interface between the pile and its surrounding soil. 3.2 Comparison of the numerical results and the measurements 3.2.1 Variation of the settlement at the centerline of the embankment Fig. 3 shows the comparison of the numerical results and the measurements for the settlement of the subsoil at the centerline of the embankment. The settlement of the subsoil increased with an increase of the embankment height. After the construction of the embankment, a few subsequent settlements occurred due to the consolidation of the foundation soil. The results of the settlement calculated using the FEM were similar to the measurements, while they showed a difference after unloading the preloading. The settlement of the subsoil decreased by approximately 4 mm for the FEM, while the change was almost negligible for the measurement. The final settlement of the subsoil derived from the FEM and the measurement was in good agreement, with a slight difference of 8.2%. The settlement of the subsoil at the centerline of the embankment can therefore be evaluated using the FEM with reasonable accuracy. Time/d 600 Figure 3. Development of the settlement of subsoil at the centerline of the embankment. 3.2.2 Settlement of the subsoil along the width of the embankment The results of the settlement of the subsoil along the width of the embankment from the FEM and the measurement are presented in Fig. 4. The maximum settlement of the subsoil occurred at the centerline of the embankment, and decreased towards the right toe of the embankment. As is clear in Fig. 4, the FEM and the measurement showed a similar profile for the settlement distribution (the curve likely to be parabolic). The settlement of subsoil, from both the FEM and the measurements in the field near Distance from the centerline/m LS 3,6 5,4 7,2 9 10,8 12,6 14,4 1С J Figure 4. Settlement of subsoil along the width of the embankment. the centerline of embankment, was in good agreement. For the measurements, the settlement at the right toe of the embankment is approximately 10.8% of that at the centerline of the soil surface, while for the FEM, the corresponding percentage is approximately 30.2%. It was concluded that the numerical method overestimated the settlement on the side of the embankment. 3.2.3 Lateral deformation of the subsoil Fig. 5 shows a comparison between the numerical results and the measurement in terms of the lateral deformation near the right toe of the embankment at the end of the construction. As shown in Fig. 5, the measured maximum lateral deformation was 7 mm, and this occurred at the surface of the subsoil, which showed a similar shape to the lateral deformation profile from the FEM. However, the magnitude of the lateral deformation was overestimated, with a difference of 17.6%, which may be due to the use of isotropic homogeneous soil models and soil parameters in the FEM to simulate what were aniso-tropic and nonhomogeneous soils in the field. Since the piles were founded on weathered rock layers rather than the fully rigid layer, the value of the lateral deformation at the bottom of the piles from the FEM was not equal to zero, which was also different to that of the measurement. Lateral derformation/mm 3.2.4 The stress-concentration ratio (SCR) The degree of stress concentration from the soil to the piles is typically evaluated using the stress-concentration ratio (SCR), which is defined as the ratio of the stress on the pile caps to that on the subsoil, and is given by: s p SCR = (3) ss where Op is the vertical stress acting on the pile caps, and os is the vertical stress acting on the subsoil. The higher the concentration ratio, the more the stress is transferred onto the piles. As is clear in Fig. 6, the value of the SCR was almost linearly increasing with an increase of the embankment height, due to the arching in the embankment, which meant that the vertical stress transferred from the subsoil to the pile caps. The value of the SCR fluctuated after the construction of the embankment for the measurement, which may be due to the consolidation of the subsoil and the influence of the rainy season. The SCR from the FEM reached a stable value after the construction of the embankment, while it decreased after unloading the preloading. It is also clear in Fig. 6 that the SCR was underestimated by the FEM, with a difference of approximately 20%. Time /d 600 Figure 5. Lateral deformations of the subsoil at the right toe of the embankment. Figure 6. Development of the stress-concentration ratios (SCR). 3.2.5 Axial force of the pile Fig. 7 shows the distribution of the axial load along the length of the pile. As is clear from the figure, the axial load was firstly increased, while it decreased after a certain depth, and the reduction ratio increased when near the bottom of the pile. This may due to the negative frictional stress existing along the pile. The results of the FEM showed a similar trend to the measurement, while they were underestimated at the top of the pile. This was consistent with the comparison of the SCR, in which the SCR from the FEM was underestimated. Less stress acting on the pile cap would result in a smaller axial force in the pile. Axial force of pile/kN 600 4. ANALYTICAL METHODS Figure 7. Axial force of the pile. 3.2.6 Frictional stress along the pile The frictional stress was generated due to the differential settlement between the piles and the surrounding soil. The negative frictional stress will be generated when the settlement of the piles was larger than that of the subsoil. In contrast, a positive frictional stress occurred. Hence, a neutral point will be generated when the settlement of the pile equals that of the subsoil, which was approximately 6.8 m depth of the pile in Fig. 8. Since the frictional stress was derived from the axial force of the pile, the comparison of the FEM and the measurement was therefore similar to that of the axial force of the pile. The negative frictional stress was underestimated by the FEM, while it was overestimated for the positive frictional stress. Figure 8. Frictional stress along the pile. As illustrated above, the results of the FEM and the measurement showed a closer agreement. The FEM gave almost the same results to that of the field measurements in terms of the settlement, but with a slight difference. The stress-concentration ratio (SCR), the axial force and the frictional stress were somehow underestimated by the FEM. The settlement of the subsoil from the FEM was therefore a reliable parameter to assess the performance of the pile embankment. The stress-concentration ratio (SCR) predicted by the modified Terzaghi method (Russell and Pierpoint, 1997) [4] will be presented in the paper, together with the settlement of the subsoil at the centerline of the embankment and the tensile force of the geosynthetic reinforcement by Zhuang et al. (2014) [10]. 4.1 Stress-concentration ratio Terzaghi (1943) [3] proposed an analytical solution based on the trapdoor experiment. Two vertical shear planes were assumed to occur in the soil as the trapdoor descended and the shear stress applied to the soil mass, directly above the trapdoor, the shear stress decreased the soil pressure on the trapdoor and increased the pressure on the adjoining stationary soil mass. It was reported by Russell and Pierpoint (1997) [4] that for the three-dimensional condition, the stress-reduction ratio (SRR) could be expressed as follows: -4haK tan ( s2-a2) SRR = - 2 2 s —a 4haK tan 1 — e (4) Based on the equilibrium equation for the vertical force at the bottom of the embankment, the SCR can be derived as follows: SCR = s2 — SRR (s2 — a2) SRRa 2 (5) where K is the empirical constant, and can be taken as 0.7, according to Terzaghi (1943) [3] and Chen et al.(2008) [16]; h is the height of the embankment; s is the pile spacing between piles; a is the width of the pile cap; and ф is the effective friction angle of the embankment. 4.2 Settlement of the subsoil at the centerline of the embankment Zhuang et al. (2014) [10] proposed an analytical model to evaluate the maximum settlement at the surface of the subsoil and the maximum strain of the geogrid. For the reinforced, piled embankment, taking the unit body immediately below the reinforcement between the pile caps, the vertical equilibrium required that: 64k„ , d (6) g_ d3 3l4 "E0l--sG = 0 hs where kg is the tensile stiffness of the reinforcement; l is the span of the tensile reinforcement, and is taken as (s-a) (1+V2)/2; S is the maximum settlement at the surface of the subsoil; hs is the thickness of subsoil; E0 is the one-dimensional stiffness, which can be derived from the Young's modulus and the Poisson's ratio; and Oq is the vertical stress at the base of the embankment, including the action of arching. The maximum settlement of the subsoil can be obtained by solving Equation 6. 4.3 Maximum tensile force in the reinforcement The maximum tensile force ( T) in the reinforcement can be expressed in terms of the maximum settlement at the surface of the subsoil, as follows: T = kg (7) g 3l2 4.4 Comparison of the results from the analytical method, the FEM and the measurements Fig. 9 shows the results of the comparison in terms of the SCR between the analytical method, the FEM and the measurements. As is clear from Fig. 9, the SCR predicted by the modified Terzaghi method increased with an increase in the embankment height h, which implied that more vertical stress was transferred onto the pile cap with an increase in the embankment height. The SCR from the modified Terzaghi method was overestimated, with a value of approximately 25%, while it was underestimated by the FEM, with a value of 20%. It showed that the modified Terzaghi method was more conservative in terms of the SCR than that of the FEM. It was also shown that the SCR was not a reliable parameter when assessing the performance of the embankment. 14 12 10 g8 Co б 4 2 -Modified Terzaghi -FEM <3D Measured 1 2 3 4 5 6 7 H / (s-a) Figure 9. Comparison of SCR with the analytical method, the FEM and the measurement. The results of the comparison for the maximum settlement of subsoil (S) are presented in Fig. 10. The maximum settlement of the subsoil increased with an increase in the embankment height, and was slightly overestimated, with a value of 5%, while the results from the FEM were underestimated, with the value of 8.2%. Both the 40 36 32 - §28 i 24 -20 -16 -Zhuang et al. (2014) -FEM t T,-' -4- Measured o 4 h/m Figure 10. Comparison of the maximum settlement of subsoil with the analytical method, the FEM and the measurement. analytical method of Zhuang et al. (2014) [10] and the FEM showed reasonable agreement with the measurement. The maximum settlement of the subsoil was therefore a reliable parameter to evaluate the performance of the piled embankment with reasonable accuracy. As shown in Fig. 11, a comparison of the maximum tensile force in the reinforcement was presented. The maximum tensile force in the reinforcement also increased with an increase of the embankment height, and was slightly overestimated, with a value of less than 7%. The results predicted by the FEM were also overestimated; however, with a value of approximately 44%. The tensile force can therefore be assessed using the analytical method presented in Zhuang et al. (2014) [10] with reasonable accuracy, while the results of the FEM were much more conservative. Figure 11. Comparison of the maximum tensile force in the reinforcement with the analytical method, the FEM and the measurement. 5 CONCLUSIONS_ A full-scale field case of the Fengyang high-speed railway was presented in this paper, in which the embank- ment was constructed over soft soil, and was improved by the combined use of the CFG pile-net composited foundation and preloading. Numerous devices were used to monitor the evaluation of the embankment. Numerical simulations and analytical methods were used to predict the behavior of the embankment. The results of the numerical simulation were compared with the measurements, in terms of the development for the settlement of subsoil at the centerline of the embankment, the distribution of the settlement of subsoil along the width of the embankment, the lateral deformation of the subsoil, the stress/concentration ratio (SCR), the axial force of the pile and the frictional stress along the pile. It showed that the results of the FEM were in good agreement with the measurements. It was also shown that the SCR was underestimated by the FEM, as well as the axial force and the frictional stress of the pile. The comparison of the analytical methods with the FEM results and the measurements also showed that the settlement of the subsoil was a reliable parameter to assess the performance of the embankment with reasonable accuracy. The SCR predicted by the modified Terzaghi method was overestimated, with a difference of 25%, while it was underestimated by the FEM, with a difference of 20%. It was also shown that the tensile force in the reinforcement can be effectively assessed by Zhuang et al. (2014) [10], while it was overestimated by the FEM, with the amount being conservative. Acknowledgements The financial support of the National Natural Science Foundation of China (Grant No. 51478166 and 51108155), the Fundamental Research Funds for the Central Universities (Grant No. 2015B06014 and 2015B17814) and the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, is acknowledged. REFERENCES_ [1] Magnan, J.P., 1994. Methods to reduce the settlement of embankments on soft clay: A review. Vertical and Horizontal Deformations of Foundations and Embankments. ASCE Geotechnical Special Publication 1, 40, 77-91. [2] Gangakhedkar R. 2004. Geosynthetic reinforced pile supported embankments. Master thesis, University of Florida. [3] Terzahgi, K. 1943. Theoretical Soil Mechanics. John Wiley and Sons, New York. [4] Russell, D., Pierpoint, N. 1997. An assessment of design methods for piled embankments. Ground Engineering 30, 10, 39-44. [5] Hewlett, W.J., Randolph, M.F. 1988. Analysis of piled embankments. Ground Engineering 21, 3, 12-18. [6] BS8006-1, 2010. Code of practice for strengthened/ reinforced soils and other fills. British Standards Institution, UK. [7] EBGEO, 2011. Recommendations for Design and Analysis of Earth Structures using Geosynthetic Reinforcements EBGEO, 2011. ISBN 978-3-43302983-1 and digital in English ISBN 978-3-43360093-1. [8] Kempfert, H.G., Göbel, C., Alexiew, D., Heitz, C. 2004. German recommendations for reinforced embankments on pile-similar elements. In Euro-Geo3-Third European Geosynthetics Conference, Geotechnical Engineering with Geosynthetics, pp. 279-284. [9] Van Eekelen, S.J.M., Bezuijen, A. and Van Tol, A.F. 2011. Analysis and modification of the British Standard BS8006 for the design of piled embankments. Geotextiles and Geomembranes 29, 3, 345-359. [10] Zhuang,Y., Wang, K.Y., Liu, H.L. 2014. A simplified model to analyze the reinforced piled embankment. Geotextiles and Geomembranes 42, 2, 154-165. [11] Han, J., Gabr, M. A. 2002. Numerical analysis of geosynthetic reinforced and pile-supported earth platforms over soft soil. Journal of Geotechnical and Geoenvironmental Engineering 128, 1, 44-53. [12] Le Hello, B., Villard, P. 2009. Embankments reinforced by piles and geosynthetics numerical and experimental studies with the transfer of load on the soil embankment. Engineering Geology 106, 1-2, 78-91. [13] Zeng, J.C. 2010. Experiment study and theoretical analysis of CFG pile composite foundation in highspeed railway. PhD thesis, Southeast University, China (In Chinese) [14] Liu, H. L., Charles, W. W., Fei, K. 2007. Performance of a Geogrid-Reinforced and Pile-Supported Highway Embankment over Soft Clay: Case Study. ASCE, Journal of Geotechnical and Geoenvironmental Engineering 133, 1483-1493. [15] Potyondy J.G. 1961. Skin friction between various soils and construction materials. Geotechnique 11, 4, 339-353. [16] Chen, Y.M., Cao, W.P., Chen, R.P. 2008. An experimental investigation of soil arching within basal reinforced and unreinforced piled embankments. Geotextiles and Geomembranes 26, 2, 164-174. SUHA FRAKCIJA NEZASIĆENIH TAL Izvleček V članku je predlagana enačba za račun strižne trdnosti nezasičenih tleh. Ta enačba je definirana kot enakovredna napetost in je razširitev Murrayeve enačbe. Ta pristop se uporablja za splošni primer dvo-modalno strukturiranih tleh, ki kažejo makrostrukturo in mikrostrukturo. Teoretični razvoj obravnava obstoj suhe frakcije poleg zasičene in nezasičene frakcije tal. Te različne frakcije so vključene v porozni model, ki omogoča ovrednotenje parametrov enačbe primerjalne napetosti. Nadalje je izvedena primerjava med teoretičnimi in eksperimentalnimi rezultati. Slednje kaže, da se predlagana enačba lahko uporabi za oceno strižne trdnosti nezasičenih tal. Julio César Leal Vaca (vodilni avtor) Universidad Autónoma de Querétaro, México, Engeenering Department 77 Juarez Ave., Downtown, C.P. 36000 Guanajuato, Gto., Mehika E-pošta: jcesarlealv@hotmail.com Gustavo Gallegos Fonseca Universidad Autònoma de Querétaro, México, Engeenering Department, 501 Romualdo del Campo, Rafael Curiel, C.P. 79060 Ciudad Valles, San Luis Potosi, Mehika E-pošta: gfonseca@uaslp.mx Eduardo Rojas Gonzalez Universidad Autònoma de Querétaro, México, Engeenering Department Centro Universitario Cerro de las Campanas, C.P. 76010, Querétaro, Qro. Mehika E-pošta: erg@uaq.mx Ključne besede nezasićena zemljina, strižna trdnost, ekvivalentna napetost, krivulja zadrževanja vode THE DRY FRACTION OF UNSATURATED SOILS Keywords unsaturated soil, shear strength, equivalent stress, water retention curve Abstract An equation to account for the shear strength of unsaturated soils is proposed in this paper. This equation is defined as the equivalent stress, and is an extension of Murray's equation. This approach applies to the general case of bi-modal structured soils showing a macrostruc-ture and microstructure. The theoretical development considers the existence of a dry fraction in addition to the saturated and unsaturated fractions of the soil. These different fractions are included in a porous model, which allows an evaluation of the parameters of the equivalent stress equation. Finally, the paper includes a comparison between theoretical and experimental results. The comparison shows that the proposed equation can be used to estimate the shear strength of unsaturated soils. 1 INTRODUCTION The volumetric behavior and shear strength of saturated soils can be clearly explained using the principle of effective stress, as established by Terzaghi [1]. For the case of unsaturated soils, a general effective stress equation has not been established yet. Bishop's equation [2] has been used as an effective stress equation for unsaturated soils. However, there is still no consensus on how to evaluate or determine the parameter x, which varies between 0 and 1, for the dry and the saturated conditions, respectively. Other effective stress equations have been proposed by different researchers, for example: Croney et al. [3], Skempton [4], Aitchison [5, 6], Jennings [7], and Richards [8]. These equations hold to the idea of effective stress and include empirical factors similar to Bishop's parameter x. For this reason, the behavior of unsaturated soils has been studied on the basis of the independent stress variables proposed by Burland [9]. In general, the net stress and the suction have been used as the two independent stressstate variables to simulate the behavior of these materials. Recently, Bishop's equation has regained importance, since new developments in the modeling of porous structures make it possible to determine the parameter X. Also, the hydro-mechanical coupling observed in unsaturated materials can be easily introduced in effective stress formulations. To that purpose, the influence of the volumetric deformation on parameter x should be included. Various equations have been proposed to determine Bishop's parameter x, for example: Öberg and Sällfours [10], Vanapalli et al. [11], Khalili and Khabbaz [12] and Rojas [13]. These equations include the degree of saturation and/or the suction as variables. Rojas [13] proposes the following equation. Julio César Leal Vaca (corresponding author) Universidad Autónoma de Querétaro, México, Engeenering Department 77 Juarez Ave., Downtown, C.P. 36000 Guanajuato, Gto., Mexico E-mail: jcesarlealv@hotmail.com Gustavo Gallegos Fonseca Universidad Autónoma de Querétaro, México, Engeenering Department, 501 Romualdo del Campo, Rafael Curiel, C.P. 79060 Ciudad Valles, San Luis Potosi, Mexico E-mail: gfonseca@uaslp.mx Eduardo Rojas Gonzalez Universidad Autónoma de Querétaro, México, Engeenering Department Centro Universitario Cerro de las Campanas, C.P. 76010, Queré-taro, Qro. Mexico E-mail: erg@uaq.mx X = fs + (1 - fs ) sw (1) where fs and SW represent the saturated fraction and the degree of saturation of the unsaturated fraction, respectively. One way to relate the suction and the degree of saturation is by means of the soil-water characteristics curves (SWCCs). There are several proposals for a general equation for the SWCCs, most of them containing empirical fitting parameters. Some of these equations have been proposed by Gardner [14], Brooks and Corey [15], Farrell and Larson [16], Van Genuchten [17], McKee and Bumb [18], Seber and Will [19] and Fredlund and Xing [20]. There are also relations oriented to determine the SWCCs from the granulometric distributions: Arya, L.M. and Paris, J. F. [21], Huang, M. [22]. Other researchers use the properties of mass and volume, such as Fredlund, M. D. [23, 24]. Sillers [25] has applied statistical methods, while Johari et al. [26] have used numerical techniques such as genetic algorithms to predict the SWCCs. Another way to obtain both the effective stress as well as the SWCC is by means of a porous-solid model. This type of model attempts to reproduce the porous structure of soils and is able to simulate the distribution of water in the pores of the soil depending on the value suction. This property makes it possible to determine Bishop's parameter x as a function of the saturated and unsaturated fractions of the soil fs, fu); for example, the the case of equation (1), in which fs + fu = 1. The saturated fraction represents all the zones where solids are completely surrounded by saturated pores, while the unsaturated fraction is formed by solids surrounded by a combination of saturated and unsaturated pores. However, in this paper it is demonstrated that a dry fraction may eventually appear at high suctions, affecting the value of the parameter x. The dry fraction is formed by solids completely surrounded by dry pores. In this paper, an equation of equivalent stress that includes the dry fraction is proposed. 2 EQUIVALENT STRESS EQUATION WITH THE DRY FRACTION_ An unsaturated soil is formed by solid particles where the voids are occupied by two other phases: air and water. The structure of the soils depends on the shape, size and stresses applied to the solid particles. The mechanical behavior of these materials is influenced by the presence of water and its interaction with air. The water contained in most soils is present in different states: as a solid in the inner layers of the so-called adsorption layer. When the molecules of water move away from the surface of the fine particles of soil, it becomes less viscous. Farther away, the water molecules show a normal viscosity where the capillary and gravitational forces are predominant. The capillary phenomenon is produced by the meniscus of water formed mainly at the contact between the solid particles. This capillary effect produces additional contact stresses between the solid particles. These additional contact stresses need to be evaluated in order to propose an effective stress equation for unsaturated materials. One way to obtain these additional contact stresses was proposed by Murray [27]. Murray applied the concept of enthalpy to determine the coupling stress between the solid particles p'c as a function of the total stress p, the air pressure ua, the value of the suction s, the void ratio e and the degree of saturation Sw, in the form f! + eSw p'c = p — ua + s 1+e (2) It is clear that equation (2) keeps the structure of equation (3) proposed by Bishop, in which X = (1 + eSw )/(1 + e). °'= s — ua + SX (3) where a' is the effective stress, a - ua is the net stress, s is the suction and x is Bishop's parameter. Following Murray's proposal, the strength of the material can be represented with the following relationship q=Mt --1 + Л (4) where Mt represents the slope of the line for the coordinated system p'c/s — q / s , and Л is the coordinate to the origin of such a straight line. Notice that equation (4) becomes undetermined for the saturated case (i.e., the suction is zero). In that sense Murray's model is unable to consider the effect of the saturated fraction of the soil. These equations were obtained by considering a monomodal structure of soils. However, according to Srid-haran et al. [28], most soils show a bimodal structure: consisting of a macrostructure and a microstructure. The first one is related to larger particles and inter-particle pores, also called macropores. The second one refers to smaller particles (for example, packets of clay) and intra-particle pores or micropores. By following Murray's procedure, but considering a bi-modal structured soil, Rojas [13] obtained an alternative expression for the equivalent stress as represented by Equation (5). This equation also follows the structure of Bishop's equation (3), where x is given by equation (1). pc= p - ua + s [ fs + SUW (1-fs )] (5) where the term (1 - fs) represents the unsaturated fraction of the soil. These unsaturated and saturated fractions naturally appear in a soil sample during a drying or wetting process. For example, at the beginning of the drying process for an initially saturated soil, all the pores are saturated and thus only the saturated fraction exists. As the suction increases, some of the largest pores start to dry and therefore some solids are surrounded by a combination of saturated and dry pores and an unsaturated fraction emerges. However, when a soil is subjected to high suctions, most of the macrostructure dries, and therefore, large particles may be completely surrounded by dry pores, and in that sense a dry fraction appears. This paper proposes an extension of Rojas's equation [13] to include the dry fraction of the soil following the procedure established by Murray to obtain a coupling stress. Fig. 1 represents the elemental volume of a bi-modal structured unsaturated soil subjected to a total stress p. The enthalpy Ht of the sample is the product pV, where У represents the volume of the elemental volume. 1ЩЩЩЩЩЩ11 Figure 1. Soil fractions under pressure. The total enthalpy of the sample is represented by the sum of the enthalpies of each one of its phases. Ht = Hu, + H„ + Hs (6) where Ht is the total enthalpy, Hw is the enthalpy of the water phase, Ha is the enthalpy of the air phase, and Hs is the enthalpy of the solid phase. Following the proposal of Murray [27], and keeping the products pV = Ht, Hw = uwVw, Ha = uaVa, Hs = usVs + p'cV, it is possible to establish the following expression: pV = uwVw + uaVa + usVs + p'cV (7) where p is the total stress, V is the total volume, uw is the water pressure, Vw is the volume of water, ua is the air pressure, Va is the volume of air, us is the stress in the solid particles, Vs is the volume of solids, and p'c is the coupled stress. The coupling stress p'c can be seen as a component of the enthalpy by unit of total volume. This stress links the independent stress variables and allows defining in a clearer way the behavior of unsaturated soils. The volume of the sample can be divided into three parts: the saturated fraction Vs, the unsaturated fraction Vu and the dry fraction Vd in the form: V = Vs + Vu + Va (8) Dividing equation (8) by V the following relationship is obtained: 1 = fs + fu + fd (9) where fs, fu and fd represent the saturated, the unsaturated and the dry fractions, respectively. Considering the unsaturated and saturated fractions according to figure 1, and using particle i as reference, equation (7) can be rewritten as follows: uwVw uwVL u"V"w uaV" uaVa Vr p = _w_w + w sw + w sw +_JLJa + _JLJL +u_SL + p'c (10) V V V V V V where V = V!w + Vl + VU + Vsc is the solids' volume, s sw sw sa sc Vsw is the solids' volume of the saturated fraction in contact with the water, VUw is the solids' volume of the unsaturated fraction in contact with the water, VUa is the solids' volume of unsaturated fraction in contact with the air, Vsc is the contact volume among the solid particles,Vw = V^ + VU , where V^ is the volume of water in the pores of the saturated fraction, and Vw is the volume of water in the pores of the unsaturated fraction, and V5C / V » 0 , according to Skempton [29]. This consideration was explained by Skempton using the area-ratio concept defined as a = Asc / A , where Asc represents the area of contact between the particles and A is the total area of the cross-section, figure 1. The total area A results from the addition of saturated As, unsaturated Au and dry Ad fractions, that is to say A = As + Au + Ad. In this area, a total normal stress p, is applied. Parameter a is so small that it can be neglected for a unit thickness Vsc / V, meaning that the contact area between the solid particles is being neglected. By grouping the terms with ua and uw in equation (10) the following is obtained p'= p - V(vw + vw + vw + vw)-V( + vu (11) where V = Vw + V^ + VI + VU + VU + Vau, then the value VUa + VU can be replaced by V-( + VU + V^ + Vsl) resulting in p'c = p - Ua +(ua - Uw) Vs + Vu + Vs + Vu w 1 w 1 sw 1 sw V (12) The volume of the saturated fraction is Vs = Vw + V^w , which when substituted into equation (12) and defining s Vs fs = results in V Vu + Vu fs +- w sw p'c = p - ua + s V (13) Reducing the second term in the parentheses of equation (13) yields the following expression. Vu + Vu w ' sw V Vu + ru Vu . w 1 sw s V suVUU + V (14) Vu Vu Vu where Suw =, ru = » bw. and - Vu ' sw vu VU Vu + Vu ,,Vu w sw = Su-v- V u w V Vu u v" u VU -si—=SU-?-w V w V V" u V" u V" u -SU— = sU—+s" w V w V w Vu -VU nU Vv I nU r nU Vv nU r nU f-, r r \ = SwV + Swfu - SwV = Swfu = Sw (1-Js - Jd ) (15) Finally, the equation for the equivalent stress that includes the dry fraction is p'c = p - Ua + s ё fs + S"w (-fs - fd ) (16) In equation (16) the value of Bishop's parameter can be easily identified as X = fs + Sw (1 - fs - fd) .For the case of a fully saturated soil, s is substituted by ua - uw, fs = 1 and fd = 0 or X = 1, and then equation (14) becomes Terzaghi's equation p'c = p - Uw (17) For the case of a dry soil, x = 0 and equation (16) transforms into p'c = p - ua (18) 3 EXPERIMENTAL PROGRAM_ In order to experimentally verify the proposed equation of equivalent stresses, a laboratory testing program was developed. It included the following steps: a) a mixture of 79% sand, 21% silt was made, b) the material index properties were determined, c) its grain size distribution was determined, d) a methodology to prepare compacted soil samples was established, e) using the filter paper technique, water-retention curves for wetting and drying were determined f) the porosimetry of the sample was obtained by fitting the numerical and experi- mental retention curves in drying and wetting, and g) triaxial tests were conducted in a controlled suction triaxial cell on samples following drying and wetting paths. Each one of these steps is explained below. a) The sand and silt were obtained from an extraction bank located in the city of Valle de Santiago, Guanajuato, Mexico. The sand was passed through the sieve number 10, and washed to remove clay particles. b) This material has neither plasticity nor linear shrinkage. The relative density of the solids is 2.43 and the soil was classified as SM (silty sand). c) All the specimens were fabricated in a metallic mould by static compaction in five layers at a water content of 19.5% and a dry density of 14.889 kN/m3. Under these conditions, the void ratio of the samples was 0.54 and their degree of saturation 87%. In order to fabricate the sample, each layer was weighed on a balance with 0.01-g precision and placed inside the mould. In order to avoid planes of contact, the surface of the previous layer was carefully scarified before placing the material for the next layer. d) The granulometry of the material was determined using the dry and wet methods. e) The size distribution of the macropores, sites and bonds was obtained by adjusting the water-retention curves during drying and wetting, taking into account the following: "it is well known that the drying curve is mainly dependent on the size of bonds as they unsaturate at higher suctions than sites. On the contrary, the wetting curve depends mainly on the size of sites and macropores as they saturate at smaller suctions than bonds. Therefore, it is proposed to use of the boundary SWCC at drying in order to define the size distribution of bonds. When the porosimetry of the material is not available at all, both retention curves have to be used to define the size distribution of bonds sites and macropores. These distributions are obtained from an iterative procedure where an initially proposed size distribution for each element is consecutively modified, until the porous model reproduces with sufficient accuracy, both the wetting and the drying soil-water characteristic curves", (Rojas et al. [13], p. 198). f) The only method to obtain directly the retention curve is with the membrane apparatus. However, its suction range only reaches 10 MPa. On the other hand, the filter-paper method covers the entire range of suctions and can be easily applied for wetting and drying paths. However, mainly due to the heterogeneity of the papers during the fabrication process, it is recommended to verify the calibration every time a new batch of paper is used. That is the reason why the Sleicher and Schuell No. 589 paper batches were first calibrated with potassium chloride solutions at various concentrations. The whole process was conducted in a controlled-temperature room where the filter paper was exposed to a saline atmosphere for a period of two weeks. Once the filter paper was calibrated, cubic samples with a side of 5 cm were made. The drying curve was obtained by introducing the compacted samples into the oven for time ranging from 1 min to 24 hours, then introduced in a sealed container along with two pieces of filter paper to measure both the total and matrix suctions. The samples were left for two weeks in a controlled--temperature room to reach equilibrium. In the case of the wetting curve the samples were initially oven dried for 24 hours and subsequently moistened using a fine spray to reach a certain degree of saturation. Once the sample reached the required humidity, the homogenization process was carried out for two weeks in an identical manner as for the drying path. g) Triaxial tests were performed following the wetting and drying paths. The suction was controlled using the vapor-circulation technique using a peristaltic pump with an 11-ml/min flow rate, as recommended by Cunningham [30]. The soil samples were placed in the triaxial chamber, and their initial moisture content was slightly higher or lower than the moisture of the circulating vapor, so as to follow a drying or wetting path, respectively. In order to define the equilibrium time for specimens, some tests were performed previously, where the samples were mounted in the triaxial chamber and weighed at regular intervals until equilibrium was achieved. In general, it was found that three or four days were sufficient to reach equilibrium. In order to ensure constant suction tests, all the tests were performed at a strain rate of 0.001 mm/min. The confining stress for all the samples was 150 kPa. Once the triaxial tests were completed, samples were cut in sections to verify the value of the suction. Triaxial tests were performed on soil samples with water contents ranging from saturated to practically dry. For the case of the saturated soil samples, confining stresses of 50, 100 and 150 kPa were applied. 4 COMPUTATIONAL MODEL_ The solid-porous model used herein is similar to that proposed by Rojas et al. [31], in which four types of elements were considered: the macropores, the sites, the bonds and the solids. This model is built on a regular network and includes the following elements: cavities, bonds, and solids. are responsible for most of the volumetric behavior of the soil (Simms and Yanful, [32]). The mesopores or sites are pores that do not change their size during loading. All the cavities are interconnected by the bonds or throats. Finally, the solids are considered incompressible. The model is built on a regular network, where the cavities are placed at the nodes of the network, and the connectors are the bonds. These elements are represented in Fig. 2. Macropare Bond Figure 2. Network elements: macropores, sites, bonds and solids (from Rojas et al. [31], p. 196). The size of the network determines the number of nodes and connectors. This defines the number of cavities, bonds and solids. The number of elements of each size is defined by the porosimetry and the grain size distribution of the material. Once the number of elements of each size is defined, sites and bonds are randomly placed on the network. In order to avoid the superposition of bonds converging on one node, a constructive principle needs to be respected. This constructive principle establishes that two bonds converging to a site at 90° should comply with expression (19). 4 b1 ' rb2 £ rst (19) where rb1 and rb2 are the radius of the converging bonds and rst is the radius of the cavity. To totally guarantee this principle within the network, this procedure must be followed: "at each node, the number of violations to the principle is determined. If the number of violations at certain site is different from zero, then a substitution with another site (selected at random) is simulated. If the number of violations to the principle reduces, then the substitution is granted. If not, another site is selected at random. The same procedure is followed for bonds, and it continues until no violation to the construction principle subsists within the network", (Rojas et al. [31], p. 197). Cavities contain most of the void volume of the soils. These cavities are subdivided into macropores and meso-pores. The macropores are the largest pores in the soil and Once the sites and bonds are located, macropores are placed by substituting the required number of sites. Finally, solids are also placed at random in the spaces between the pores. The process of drying and wetting the pores is regulated by the Young-Laplace equation (20), Defay and Prigogine [33]. This equation is a nonlinear partial differential equation that describes the capillary pressure difference sustained across the interface between two static fluids, such as water and air, due to the phenomenon of surface tension or wall tension, although usage on the latter is only applicable if we assume that the wall is very thin. The Young-Laplace equation for a spherical meniscus can be written as: 2T cos a (20) u - u, = - where a represents the contact angle between the water and the solid particles, Rc is the critical radius and corresponds to the maximum pore size that remains saturated at a certain suction. In addition, to comply with equation (20) a pore must also comply with the continuity principle in order to saturate or dry. This continuity principle establishes that a pore is able to dry or saturate only if it is connected to a boundary of the network where the bulk of gas or water is present, respectively, following a continuous path. Therefore, during a drying process, pores will dry when the size of the pore is larger or equal to Rc and one of its connected bonds is already dry and connected to the bulk of gas. The drying process starts when all the pores are saturated and the suction is zero. An increase in the suction results in the drainage of the largest pores located at the boundaries of the network. The process continues only when the bonds connected to these pores dry too. This means that the drying SWCC is dependent on the size distribution of the bonds. The saturation of pores will happen only if the pore size is smaller or equal to Rc and one of the elements connected to this pore is already saturated. During the wetting process it is considered that initially all the pores are dry and the suction is very high. Then the suction decreases in steps and the smallest bonds at the boundaries of the network start to saturate. The wetting process only continues when the sites connected to those bonds saturate. This means that the wetting process is controlled by the size distribution of the sites. A computer program was developed to reproduce the porous structure of the soils by fitting the numerical with the experimental SWCC. The program also determines the saturated, unsaturated and dry fractions of the soil during wetting-drying cycles. Finally, the parameter x is determined and the experimental and numerical results are compared. In general, logarithmic normal distributions are used to define the pore size distributions of soils. Logarithmic normal distributions are defined with only two parameters: the mean value and the standard deviation. Therefore, these two parameters are required for each of the elements in the network: sites, macropores, bonds and solids. Double logarithmic normal distributions are required for the case of double structured soils. 5 DISCUSSION OF THE THEORETICAL AND EXPERIMENTAL RESULTS._ Table 1 indicates the parameters required by the solid-porous model: the mean radius ( R ) and the standard deviation (S) for sites (SPX and SP2), macropores (MP) and bonds (Bx and B2). Finally, the void ratio (e) of the sample is also required. Table 1. Parameters of the solid-porous model. Parameter SP1 SP2 MP B1 B2 e R 0.05 0.4 3.0 0.0004 2.0 0.54 S 2 2 3 6 4 The fitting of the numerical with the experimental SWCC is shown in Fig. 3. The values included in table 1 define the logarithmic normal distribution for each case. These parameters are initially proposed and then adjusted until the best fit with the experimental points is achieved. The experimental points of the SWCCs were obtained using the filter-paper method. -Theoretical (drying) --Theoretical (wetting) * Experimetal (drying) ■ Experimental (wetting) o,ooood,oooi 0,001 0,01 0,1 s(MPa) 1000 Figure 3. Numerical and experimental SWCC in drying and wetting paths (data from Leal et al. [34]). The numerical pore size distribution obtained after the fitting process is shown in Fig. 4. The existence of a bimodal pore structure can be observed. 1000 10000 Radius ( um j Figure 4. Theoretical pore size distribution of the soil. In the same way, in order to obtain the number of solids of each size, the experimental grain size distribution has to be fitted with the numerical model using a single or double logarithmic normal distribution. These results are shown in Fig. 5. It is clear that the fitted curve is quite close to the experimental results. -Theoretical Eperimentäl £ I Fh v? O'- 1000 Gain 51яе (|_aii) Figure 5. Theoretical and experimental grain size distribution. - Drying-CÜJ Д Drying-L/i) - Drying-W) Wetting-[/d) Wetting-( Я) Wetting-Ш О 0.1 02 03 04 0.5 0.6 0.7 OS 05 1 Sw Figure 6. Theoretical values of the dry, unsaturated and saturated fractions for the drying and wetting paths. -----Drying -Wetting 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 Sw Figure 7. Theoretical values for the degree of saturation of the unsaturated fraction in drying and wetting (from Leal et al. [34]). With the model, the parameters fs, fu, fd, and Swu for different degrees of saturation (Sw) are defined. These parameters are shown in Fig. 6 and Fig. 7. The volume of the saturated fraction ( fs = 1 — fu — fd ) is obtained by adding the volume of solids completely surrounded by water to the volume of voids surrounding these solids. Fig. 6 shows the values of the dry fraction (fd), the unsaturated fraction (fu), and the saturated fraction (fs) with the degree of saturation (Sw). These values were obtained for the drying and wetting paths. It can be verified that the values shown in Fig. 6 comply with Equation (9) for any value of the degree of saturation for both the drying and wetting paths. Additionally, Fig. 7 shows the variation of the degree of saturation for the unsaturated fraction with respect to the degree of saturation for both paths. Bishop's parameter can now be determined using the expression c = fs + Sw (1 — fs — fd) included in Equa- tion (16). To obtain the experimental values of x, the p — q diagram is plotted with the triaxial test results of saturated samples, carried out at confining pressures of 50, 100 and 150 kPa. According to the critical state theory it is possible to draw a failure line departing from the origin of the plane p — q. If in this diagram the pairs of values (ffnet, q) of the unsaturated samples are represented, then the difference between their abscissas and the failure line represents the suction stress sx. And because the suction is known, it is possible to obtain the value of x. The experimental and numerical values for x are shown in Fig. 8, for both the drying and wetting paths. In this figure, it can be observed that both the numerical and experimental values follow a general tendency, except for one experimental point obtained during wetting. "The best results are presented in low degrees of saturation. This is because at low degrees of saturation volume changes are minor and the degree of saturation is influenced by changes in volume. One can expect some differences in the final Theoretical (drying) Experimental (drying) -Theoretical (wetting) A Experimental (wetting) 1 -i 0,9 ■ o,s ■ 0.7 ■ 0.5 -I 0.4 0.3 0.2 ■ 0.1 ■ 0 0 0.1 0.2 0,3 0.4 0.5 0,6 0.7 0,8 0,9 1 Figure 8. Variation of x with respect to the degree of saturation (data extended, from Leal et al. [34]). 900 soo "00 600 _ 500 3, 400 300 200 100 0 100 200 300 400 p' (kPa) 500 600 Figure 9. Numerical results of the shear strength using the solid-porous model. results off and therefore also of the resistance because the processes of wetting and drying of the solid-porous model are developed by invasion. That is, to move water or gas, continuity in the phases is required. This means that a site or a bond cannot be invaded if at least one adjacent element has not been invaded", (Leal et al. [34], p.401). When the value of x is multiplied by its corresponding suction, the cohesive stress is obtained. Once the values of the parameter x have been determined, it is possible to estimate the equivalent stress values using the following expressions: 0-1 = 03 + q (21) Si neto 0ì ua i neto (22) (23) In these equations, i varies from 1 to 3, a1 is the total vertical stress, a3 is the confining stress, ua is the air pressure in the pores, q is the deviator stress applied in the triaxial tests, 0' is the equivalent stress and s is the suction. All the specimens were tested up to failure according to the critical state theory. The relationship between the mean and the deviatoric stress is represented by the slope M. If the results are plotted in this plane, their alignment with the slope M can be observed, and the mean equivalent stress can be obtained using the following expression: po = (o1 + 02 + 03) / 3 (24) Following the above equations, it is possible to calculate the theoretical equivalent stress pc', and the results are shown in Fig. 9. Fig. 9 shows that the solid-porous model can be used to obtain Bishop's parameter x and simulate the strength of the unsaturated soils tested at different values of suction. Finally, the inclusion of the dry fraction for the determination of Bishop's x parameter results in a more precise description of the behavior of unsaturated soils during the wetting-drying cycles. 6 CONCLUSIONS_ A general analytical equation has been established to obtain Bishop's parameter x, that includes the saturated, the unsat-urated and the dry fractions of the soil. The comparisons of the numerical and experimental results show that the proposed equation is adequate for simulating the strength of unsaturated soils. Although more comparisons are still required, the solid-porous model proposed herein to obtain the parameter x seems to be sufficiently accurate. REFERENCES [1] Terzaghi, K. 1936. The shear resistance of saturated soils. In: Proc. 1st Int. Conf. Soil Mech. Found Eng. Cambridge, M.A. 1, pp. 54-56. [2] Bishop, A.W. 1959. The principle of effective stress. Teknisk Ukebald 106, 859-863. [3] Croney, D., Coleman, J. D., Black, W. P. M. 1958. Movement and distribution of water in soil in relation to highway design and performance. Water and its conduction in soils, highway res. Board, special report. Washington, DC. 40, pp. 226-252. [4] Skempton, A. W. 1960. Terzaghi's discovery of effective stress. In: From theory to practice in soil mechanics. John Wiley, New York. [5] Aitchinson, G. D. 1961. Relationship of moisture and effective stress functions in unsaturated soils. In: Pore pressure and suction in soil Conf., organized by British Nat. Soc. of Int. Soc. Soil Mech. Found. Eng. at Inst. Civil Eng. London, England. Butterworths, pp. 47-52. [6] Aitchinson, G. D. 1965. Soil properties, shear strength, and consolidation. In: Proc. 6th Int. Conference Soil Mech. Found. Eng. Montreal Canada. 3, pp. 318-321. [7] Jennings, J. E. 1961. A revised effective stress law for use in the prediction of the behavior of unsatu-rated soils. Pore pressure and suction in soils Conf., organized by British Nat. Soc. of Int. Soc. Soil Mech. Found. Eng. at Inst. Civil Eng. London, England. Butterworths, pp. 26-30. [8] Richards, B. G. 1966. The significance of moisture flow and equilibria in unsaturated soil in relation to design of engineering structures built shallow foundations in Australia. In: Symp. on permeability and capillary, Amer. Soc. Testing. Materials. Atlantic City, NJ. [9] Burland, J. 1964. Effective stresses in partly saturated soils. In: discussion of Some aspects of effective stress in saturated and partly saturated soils, by G. E. BlightIn, Moist and A. W. Bishop, Geotechnique 14, 65-68. [10] Oberg, A., Sallfours, G. 1997. Determination of shear strength parameters of unsaturated silts and sands based on the water retention curve. Geotech-nical Testing Journal 20, 40-48. [11] Vanapalli, S.K., Fredlund D.G., Pufahi D.E., Clifton A.W. 1996. Model for the prediction of shear strength with respect to soil suction. Canadian Geotechnica Journal 3, 379-392. [12] Khalili, N., Khabbaz, M.H. 1998. Unique relationship for the determination of the shear strength of unsaturated soils. Geotechnique 40, 5, 681-687. [13] Rojas, E. 2008. Equivalent stress equation for unsaturated soils. I: Equivalent stress. International Journal of Geomechanics ASCE 8 (5), 285-290. [14] Gardner, W.R. 1958. Some steady state solutions of the unsaturated moisture flow equation whit application to evaporation from a water-table. Soil science 85, 228-232. [15] Brooks, R.J., Corey, A. T. 1964. Hydraulic properties of porous media. Hydrol. Colo. State univ., Fort Collins. March. 27, 3. [16] Farrel, D. A., Larson W. E. 1972. Modelling the pore structure of porous media. Water Resow: Res. 3, 699-706. [17] Van Genuchten, M. T. 1980. A closed-from equation for predicting the hydraulic conductivity of unsaturated soil. Soil science of America Journal 44, 892-898. [18] McKee, C.R., Bumb, A.C. 1984. The importance of unsaturated flow parameters in designing a monitoring system for hazardous wastes and environmental emergencies. Proceedings, Hazardous Materials Control Research Institute National Conference, Houston. TX, March 1984, pp. 50-58. [19] Seber, G. A. P., Wild, C. J. 1989. Nonlinear regres- sion. John Wiley & Sons, Inc., New York, N.Y. [20] Fredlund, D., Xing A. 1994. Equations for the soil-water characteristic curve. Canadian Geotechnical Journal 31, 3, 521-532. [21] Ayra, L.M., Paris, J. F. 1981. A physico empirical model to predict the soil moisture characteristic from particle-size distribution and bulk density data. Soil. Sci. Soc. AM. J. 45, 1023-1030. [22] Huang Minghin, Fredlund, D.G., Fredlund, M.D. 2009. Estimation of SWCC from grain-size distribution curves for loess soils in China. GeoHalifax. [23] Fredlund, M.D., Wilson, G.W., Fredlund D.G. 1997. Prediction of the soil-water characteristic curve from grain-size distribution and volumemass properties. Brazilian symposium on unsatu-rated soils, Rio de Janeiro, pp 22-25. [24] Fredlund, M.D., Wilson, G.W., Fredlund, D.G. 2002. Use of the grain-size distribution for estimation. Canadian Geotechnical Journal 39, 1103-1117. [25] Siller, S.W., Fredlund D.G. 2001. Statical assessment of soil-water characteristic curve models for geotechnical engineering. Canadian Geotechnica Journal 38, 1297-1313. [26] Johari, A., Javadi, A. A., Habibagahi, G. 2011. Modelling the mechanical behaviour of unsaturated soils using a genetic algorithm-based neural network. Computers and Geotechnics 38, 2-13. [27] Murray, E. J. 2002. An equation of State for Unsatu-rated Soils. Canadian Geotechnical Journal 39, 125-140. [28] Sridharan, A., Altschaeffl, A.G., Daimond, S. 1971. Pore size distributions studies. J. Soil Mech. Fond. Div., ASCE 97(5), 771-787. [29] Skempton, A. W. 1961. Effective stress in soils, concrete and rocks-Pore pressure and suction in solis. Conference of the British National Society. Butterworths London. [30] Cunningham, M.R., Ridley A.M., Dinnen, K., Burland, J.B. 2003. The mechanical behaviour of a reconstituted unsaturated silty clay. Geotechnique 53, 183-194. [31] Rojas, E., Zepeda, A., Perez., M.L., Leal, J., Gallegos, G. 2009. A four elements porous model to estimate the strength of unsaturated soils. Geotech. Geol. Eng. Published on line (October). [32] Simms, P.H., Yanful, E.K. 2001. Measurement and estimation of pore shrinkage and pore distribution in a clayey till during soil-water characteristic curve tests. Canadian Geotechnical Journal 38, 741-754. [33] Defay, R., Prigogine, I. 1966. Surface Tension and adsorption. Longmans Green, London. [34] Leal, J.C., Gallegos, G., Rojas E. 2012. The decrease of the strength of unsaturated silty sand. Ingenieria Investigación y Tecnologia XIII, 4, 393-402. KOROZIJSKI MEHANIZMI CEMENTIRANIH ZEMLJIN ZA TRI RAZLIČNE RAZTOPINE SULFATOV Izvleček Izveden je bil niz preizkusov na blokih cementiranih zemljin izpostavljenih različnim koncentracijam raztopin H2SO4, MgSO4 in Na2SO4 , z namenom simulacije in preučitve učinkov korozije na tlačno trdnost cementiranih zemljin, ki so lahko izpostavljene onesnaženemu okolju. Rezultati preizkusov kažejo, da se stopnja korozije splošno poveča s časom korozije in koncentracije raztopine, medtem, ko tlačna trdnost pada z naraščajočo stopnjo korozijo. Stopnja korozije je največja za raztopino Na2SO4, kateri sledita raztopini MgSO4 in H2SO4. Namreč, če v raztopini obstajajo SO42- ioni, stopnja korozije za pozitivne ione sledi v padajočem redu Na+, Mg2+ in H+. Po koroziji in ionski koncentraciji je bila z rentgenskimi žarki izvedena difrakcijska fazna analiza na vzorcih cementiranih zemljin. Rezultati kažejo, da tlačna trdnost pada z večanjem koncentracije Mg2+ v raztopini MgSO4 in koncentracije Na+ v raztopini Na2SO4. Medtem, ko se trdnost veča z večanjem pH vrednosti raztopine H2SO4. Na osnovi rezultatov kemijskih analiz, se korozija raztopin H2SO4 ali MgSO4 na cementiranih zemljinah razume kot sestavljen učinek, ki vključuje kombinirane korozijske procese razpadanja in kristalizacije medtem, ko korozijo raztopine Na2SO4 sestavlja učinek raztapljanja in kristalizacije. Pengju Han Taiyuan University of Technology, College of Architecture and Civil Engineering Taiyuan, Shanxi 030024, Kitajska Chao Ren Taiyuan University of Technology, College of Architecture and Civil Engineering Taiyuan, Shanxi 030024, Kitajska Xiaohong Bai Taiyuan University of Technology, College of Architecture and Civil Engineering Taiyuan, Shanxi 030024, Kitajska Y. Frank Chen (vodilni avtor) The Pennsylvania State University, Department of Civil Engineering Middletown, PA 17057, ZDA E-pošta: yxc2@psu.edu Ključne besede cementirana zemljina, tlačna trdnost, mehanizem korozije, sulfat, onesnaževanje, raztopina CORROSION MECHANISMS FOR CEMENTED SOILS IN THREE DIFFERENT SULFATE SOLUTIONS Abstract In order to simulate and study the corrosion effects on the compressive strength of cemented soils that could be exposed in a polluted environment, a series of tests were conducted on cemented soil blocks cured with different concentrations of H2SO4, MgSO4, and Na2SO4 solutions. The test results show that the corrosion degree generally increases with the corrosion time and the solution concentration, while the compressive strength decreases with the increasing corrosion degree. The corrosion degree is highest for the Na2SO4 solution, followed by the MgSO4 and H2SO4 solutions. Namely, when the SO42- ion exists in a solution, the corrosion degree for the positive ions follows this descending order: Na+, Mg2+, and H+. X-ray diffraction (XRD) phase analyses were performed for the cemented soil samples after corrosion and ionic concentrations. The results show that the compressive strength decreases with an increase of the Mg2+ concentration in the MgSO4 solution and the Na+ concentration in the Na2SO4 solution. At the same time, the strength increases with an increase of the pH value of the H2SO4 solution. Based on the chemical analysis results, the corrosion of H2SO4 or MgSO4 solutions on cemented soils is deemed as a composite action involving the combined resolving and crystallizing corrosion processes. Furthermore, the corrosion of the Na2SO4 solution of cemented soil is a composite action consisting of dissolving and crystallizing. Pengju Han Taiyuan University of Technology, College of Architecture and Civil Engineering Taiyuan, Shanxi 030024 P. R. China Chao Ren Taiyuan University of Technology, College of Architecture and Civil Engineering Taiyuan, Shanxi 030024 P. R. China Xiaohong Bai Taiyuan University of Technology, College of Architecture and Civil Engineering Taiyuan, Shanxi 030024 P. R. China Y. Frank Chen (corresponding author) The Pennsylvania State University, Department of Civil Engineering Middletown, PA 17057, USA E-mail: yxc2@psu.edu Keywords cemented soil, compressive strength, corrosion mechanism, sulfate, pollution, solution 1 INTRODUCTION The cemented soil technique is a method of mixing the cement with in-situ soils in order to improve the soil's properties. Cement-stabilized soils may or may not work in corrosive conditions. Cemented soils would be influenced under the environment of acid rain, seawater invasion, or industrial pollution. These adverse effects could result in a structural deterioration. At a corrosive site, the strength of the soil stabilized by the cement is increased at the beginning of the stabilization, but it will be decreased due to the deterioration over time. Cemented soils are utilized in a SO42- corrosive environment when the groundwater is polluted or used under the seawater. Such a corrosive environment will inevitably corrode the cemented soil and thus change its mechanical properties. This serious issue must be considered in any practical project. Several researchers have looked at the influence of a corrosive environment on the properties of cemented soils for various aspects. Venkatarama and Jagadish [1] performed an experimental study on the influence of soil composition on the strength and durability characteristics of soil-cement blocks. Walker [2] summarized the results from a comprehensive investigation undertaken to assess the influence of soil characteristics and cement content on the physical properties of stabilized soil blocks. Pei and Yang [3] stated that the influence of Na2SO4 on cemented soils contributed to the diffusion of Na+ and SO42- and that their chemical reaction with Ca(OH)2 weakened the strength of the cemented soils due to decomposition and breaking. Shihata and Bagh-dadi [4] investigated the durability characteristics and compressive strength of cemented soils after a prolonged exposure to saline ground water. Jiao and Liu [5] analyzed the influence of the cement-soil strength in an acid environment. Kolias et al. [6] studied the effectiveness of high-calcium fly ash and cement in stabilizing fine-grained clayey soils in the laboratory. Ning et al. [7, 8] investigated the behaviors of cemented soils under various environmental conditions and concluded that, in contrast to the mechanical strength, the environmental corrosion had little effect on the fracturing process. Dong et al. [9, 10] tested the mechanical properties and electrical resistance of cement-soil polluted by H2SO4. Iyengar and Al-Tabbaa [11] presented two effective magnesia-based cements for the stabilization/ solidification (S/S) of contaminated soils. Xing et al. [12, 13] showed that Mg2+, Cl-, and SO42- caused not only a change in the microstructures of salt-rich soil-cements, but also reduced the strength of the soil-cement composite. Zandieh and Yasrobi [14] proved that two polymer materials could be used to stabilize the polluted soils on road shoulders, slopes, and military and airport pads. Heineck et al. [15] analyzed the microstructural behavior of composite mixtures of residual soils and sodic bentonites that were used as contaminant barriers. They carried out a series of microstructural analyses, including X-ray diffraction (XRD), scanning electronic microscopy (SEM), and energy-dispersive spectrometry (EDS). Voglar and Lestan [16] set up a strength model that the formulations of ordinary Portland cement (OPC), calcium aluminate cement (CAC) and pozzolanic cement (PC) and additives were used for the solidification/stabilization (S/S) of soils from a contaminated industrial brownfield. Liu et al. [17] deduced the relationship between the ion concentration and the corrosion time based on the theory of chemical dynamics, damage mechanics, and the chemical reaction formula between MgCl2 and the cement soil. Yang et al. [18, 19] considered the factors of cement content, curing age concentrations of magnesium sulfate, and pH value and concluded that the strength of the cemented soil increased with the cement content and curing age. They also analyzed the microscopic mechanism of failure. Nevertheless, studies on the influence of different sulfate solutions on the properties of cemented soils are still rather limited. To further investigate the corrosion effect and process in a sulfate corrosion environment, a series of tests, including unconfined compressive tests, were conducted on the cemented soil blocks that were cured by H2SO4, MgSO4, and Na2SO4 solutions with different concentrations. Photos of the blocks' appearances, XRD phase of the corrosion power of the cemented soils, and the concentrations of SO42-, Mg2+, H+ or Na+ in the corrosive solutions after curing the blocks were also measured. Corrosion mechanisms for different sulfate solutions on cemented soils were investigated. 2 EXPERIMENTAL PROCEDURES 2.1 Materials The main chemical composition of the cemented soils for testing is listed in Table 1. A summary of the materials and the mixing machine is provided as follows. (1) Soil: Air-dried silt soil with plasticity index (Ip) = 8.1, uniformity coefficient (Cu) = 26.67, and curative coefficient (Cc) = 1.35. (2) Cement: Ordinary Portland cement (OPC) with compressive strength = 32.5MPa after 28 days of curing, produced by a local cement company in Taiyuan, a northern city in China. (3) Proportion of cemented soil contents: soil in mass: cement: water = 100:20:50. Tap water is used for the cemented soil, and water content value in cemented soil = 29.4%. The cement content value is 11.8%. (4) Mixing machine: HJW-30 blender having a body volume = 30L and a rotational speed = 48 rpm. (5) Typical size of scaled cemented soil samples: 70.7 x 70.7 x 70.7 mm3. (6) Standard curing time for the cemented soil prior to submersion with a sulfate solution: 7 days. Table 1. Main chemical composition of cemented soils. pH Ca2+ Mg2+ SO42- Cl- CO32- OH- (-) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) 9.80 427.84 354.17 2879.10 282.32 1090.8 496.64 2.2 Preparations of H2SO4, MgSO4, and Na2SO4 solutions Based on the Chinese national standards GB 50021 (2009) [20] and GB 50046 (2008) [21], the following concentrations of H2SO4, MgSO4, and Na2SO4 sulfate solutions were considered in this study: 1.5g/L, 4.5g/L, 9.0g/L, and 18.0g/L. The corrosion results for the three solutions are summarized in Tables 2-4, in which the corrosion ratings were based on GB 50021 (2009) [20]. According to GB 50021 (2009) [20], the circumstance would be graded as B for samples fully cured by pure water. For either solution with lower concentration, the corrosion rating is consistent as weak (for 1.5g/L concentration) or medium (for 4.5g/L concentration), Tables 2-4. For higher concentration (9.0 or 18.0g/L), no consistency for the corrosion rating is noted. Table 2. Chemical ingredients of H2SO4 solution and corrosion-evaluation results. Concentration pH SO 24 Corrosion Corrosion Corrosion Corrosion (g/L) degree degree (g/L) rating (g/L) rating 1.5 5.4 weak 1.302 weak 4.5 4.3 medium 2.401 medium 9.0 2.0 strong 3.358 strong 18.0 1.6 strong 3.889 strong Table 3. Chemical ingredients of MgSO4 solution and corrosion-evaluation results. Concentration Mg2+ SO 24 Corrosion Corrosion Corrosion Corrosion (g/L) degree degree (g/L) rating (g/L) rating 1.5 2.196 weak 0.960 weak 4.5 2.293 weak 2.000 medium 9.0 2.368 weak 2.168 medium 18.0 2.538 weak 2.500 medium Table 4. Chemical ingredients of Na2SO4 solution and corrosion-evaluation results. Concentration Na+ SO 24 Corrosion Corrosion Corrosion Corrosion (g/L) degree degree (g/L) rating (g/L) rating 1.5 0.521 no 1.001 weak 4.5 1.198 no 2.300 medium 9.0 1.461 no 2.805 medium 18.0 1.670 no 3.207 strong 2.3 Testing Procedure The testing procedure consisted of the following steps. (1) Prepare the cemented soil blocks, as described in the above. (2) Cure the blocks with MgSO4, H2SO4, and Na2SO4 solutions in five different concentrations, i.e., 0, 1.5, 4.5, 9.0 and 18.0g/L. (3) Photograph each soil block at the curing times of 3, 7, 14, and 28 days to observe the change in appearance during the corrosion process. (4) Conduct the unconfined compression tests for the blocks at the same curing time as described in Step 3. For each curing time, three blocks were tested simultaneously for each prescribed concentration in a solution. The average value from the three tests is used as the block strength. (5) Measure the main ionic concentrations of the solutions immediately after the blocks are removed. (6) Perform phase analyses for the cemented soil samples in powder form after corrosion using the TD 3500 X-ray diffraction (XRD) machine, made in Denmark. The cemented soil powder samples were taken from MgSO4, H2SO4, and Na2SO4 solutions with a concentration of 18.0g/L. The XRD test was a continuous scan with a scanning angle of 200-700 and a scanning speed of 0.020/s. After the test, the corrosive powders were analyzed using JADE5.0 software to determine their chemical components [22]. 3 RESULTS AND DISCUSSION_ 3.1 Change of block appearance Figures 1-3 show the photographed appearances of the samples after 28 days of curing in the solutions. From these photos, it is clear that the sulfate solution changes the look of a cement-soil block by ways of peeling, size reduction, and cracking. Detailed observations are described as follows. (1) As shown in Figure 1, for the samples cured in the H2SO4 solution the main phenomena are peeling, size reduction, and increasing corrosion degree with the greater concentration. (2) Referring to Figure 2, for the MgSO4 solution with a lower concentration (i.e., 1.5-4.5 g/L), only crystal-like material was seen on the surface. The change in the surface size is not obvious. However, for a higher concentration (i.e., 9.0-18.0g/L), the corrosion degree becomes more serious so that both crystal-like material and peeling were observed. The corrosion degree is less than that in the H2SO4 solution. At the onset of the surface peeling, the soil sample started to crack at a concentration of 18.0g/L. (3) For samples in the more dilute Na2SO4 solution with the lower concentration (1.5-4.5g/L), similar to the MgSO4 solution only crystal-like material was observed on the surface (Figure 3). For a higher concentration (9.0-18.0g/L), the corrosion degree is more significant and the sample begins to break. For samples cured in the solution with the concentration of 18.0g/L, the damage was too severe to perform a strength test. (a) 1.5 g/L (b) 4.5 g/L (c) 9.0 g/L (d) 18.0 g/L Figure 1. Photos for the sample blocks in the H2SO4 solution after 28 days of curing. (a) 1.5 g/L (b) 4.5 g/L (c) 9.0 g/L (d) 18.0 g/L Figure 2. Photos for the sample blocks in the MgSO4 solution after 28 days of curing. (a) 1.5 g/L (b) 4.5 g/L (c) 9.0 g/L (d) 18.0 g/L Figure 3. Photos for the sample blocks in the Na2SO4 solution after 28 days of curing. Based on the above observations, the influence of the sulfate solution on the soil blocks increases with the increase of the solution concentration and the curing time. For inorganic compounds with a high solution concentration (9.0-18.0g/L), the influence degree for the three solutions follows such a descending order: Na2SO4 > MgSO4 > H2SO4. Namely, when the SO42- ion exists in a solution, the positive Na+ ion has the highest corrosion degree, followed by Mg2+ and H+. 3.2 Unconfined compressive-strength test results The compressive strength of the soil block cured in a sulfate solution (f'cu) is determined by Equation 1. f = a f J cu J a (1) where a is the modified coefficient (Table 5) reflecting the corrosion degree and fcu is the compressive strength of the block cured in pure water (i.e., when a = 1). The calculated unconfined compression strengths for the cemented soil blocks are shown in Figure 4. As demonstrated, the compressive strength decreases with the increase of the solution concentration. Generally, it increases with the corrosion time, except for the Na2SO4 solution with the high concentration (18.0g/L), where the strength is reducing to zero as the sample becomes broken. Coefficients a greater than 1 are indicated in boldface in Table 5. For the a coefficients corresponding to 28 days of curing, we observed the following findings from Table 5. Table 5. Modified coefficients (a) of cemented soils. Concentration H2SO 4 MgSO4 Na2SO4 (g/L) 3 d1 7 d 14 d 28 d 3 d 7 d 14 d 28 d 3 d 7 d 14 d 28 d 1.5 1.00 0.87 0.89 0.78 1.25 1.24 1.14 1.13 1.30 1.21 1.19 1.17 4.5 0.99 0.82 0.88 0.70 1.12 1.03 0.95 0.87 1.15 1.12 1.07 0.97 9.0 0.99 0.79 0.84 0.64 1.00 0.98 0.83 0.63 0.92 0.9 0.83 0.53 18.0 0.98 0.70 0.55 0.49 0.97 0.92 0.78 0.56 0.90 0.85 0.44 0 1 d = days (1) For the MgSO4 and Na2SO4 solutions with a 1.5g/L concentration, the coefficient a is greater than 1. This is a favorable environment as the compressive strength of the cemented soils will be increased. (2) For all the other concentrations of either solution the a coefficient is less than 1, indicating that the degree of corrosion is higher and the strength is reducing. The worst case occurs when the sample is cured in the Na2SO4 solution with a 18g/L concentration, where the strength becomes zero (i.e., a = 0) due to broken samples (Table 5). For safer, optimal, and more economical designs, it is suggested that the existing environmental condition of the cemented soils be taken into account when determining the foundation bearing capacity and settlement. 3.3 Relationship between the compressive strength (У and the ion concentration of solution (C) (1) fcu versus C of SO42- The relationship between fcu and C of the SO42- ion is shown in Figure 5. As indicated, the test data seem to be encompassed by a triangle with the slant described by fcu = -0.0008C + 4. This signifies the close relationship between the compressive strength and the SO42- concentration. 3.5 2.8 n 2.1 :> 1.4 0.7 0.0 "Og/L ■ 1.5g/L ■4.5g/L "9.0g/L ■18.Og/L 7 14 21 28 Time (days) (a) H2SO4 solution 7 14 21 Time (days) (b) MgSO4 solution ■Og/L ■1.5g/L ■4.5g/L ■9.Og/L " 18.Og/L 14 21 : Time (days) (c) Na2SO4 solution Figure 4. Relationship between the compressive strength (fcu) of the cemented soils and the corrosion time. 5,0 4,0 3,0 2,0 1,0 0,0 fcu = -0.0008C +4 0 t О юоо 2000 C (mg/L) 3000 4000 Figure 5. Relationship between fcu and C of SO42-. (2) fcu versus the concentration of a positive ion (Mg2+, H+, or Na+) The relation between fcu and the positive ion concentration or the pH values is shown in Figure 6. It is clear that the compressive strength (fcu) decreases with the increase of the Mg2+ concentration (C) of the MgSO4 solution, nonlinearly expressed by fcu = -0.136ln(C) + 2.1322, where fcu is in MPa and C is in mg/L (Figure 6a). fcu decreases with the increase of the Na+ concentration of the Na2SO4 solution linearly described by fcu = -0.0006C + 2 (Figure 6c). Again, a close relationship between the compressive strength and the ion concentration is noted. fcu increases with an increase of the pH value of the H2SO4 solution, also linearly represented by fcu = 0.0897У + 0.7372, where У is the pH value (Figure 6b). The above-derived mathematical expressions for fcu could be useful for practical designs. 3.4 Analysis of the cemented soil phases and discussion of the corrosive mechanisms After 28 days of curing, XRD tests were performed for the corrosive powder samples that were taken from the MgSO4, H2SO4, and Na2SO4 solutions with a concentration of 18.0g/L. At the conclusion of the testing, the powder chemical products (phases) were analyzed using JADE5.0 software [22]. The analysis results are presented in Figures 7-9. As seen from the figures, the resulting chemical products are very different among the different sulfate solutions and pure water. Therefore, the chemical mechanism can be determined based on the X-ray diffraction (XRD) results and the chemical formulas. (1) Corrosion mechanism for the cemented soil in the H2SO4 solution The main chemical reaction for the cemented soil in the H2SO4 solution is expressed by Equations (2)-(4). (2) Ca(OH)2 + H2SO4 = CaSO4-2H2O 3CaO-2SiO2-3H2O + 3H2SO4 = 3[CaSO4-2H2O] + 2SiO2+ 6H2O (3) 3CaO-Al2O3-3H2O + 3H2SO4 = 3[CaSO4-2H2O] + Al2O3 + 6H2O (4) As indicated in the chemical formulas (2)-(4), H2SO4 reacts with Ca(OH)2 3CaO-2SiO2-3H2O (C-S-H) and 3CaO-Al2O3-3H2O(C-A-H), where H+ participates actively and causes the cemented soil to form an unsteady structure. Therefore, the corrosion of the cemented soils in the H2SO4 solution is deemed as being a "resolving corrosion". 2>5 2,0 Q_ 1,5 2 1,0 0,5 0,0 fai= -0.1361n(C) +2.1322 О 0 600 2400 1200 1800 C (mg/L) (a) Compressive strength (fcu) versus Mg2+ concentration (C). fcu = 0.0897V + 0.7372 2,5 - 2,0 - 1,5 - 2 4-3 1,0 - 0,5 ■ 0,0 ■ 2 4 6 8 10 12 pH Value (b) Compressive strength (fcu) versus pH value (V). 2,5 -i 2,0 - 1,5 - 1,0 ■ Ч.У 0,5 ■ 0,0 ■ О О О fcu = -0.0006C +2.13 200 600 1000 1400 1800 2200 C (mg/L) (c) Compressive strength (fcu) versus Na+ concentration (C). Figure 6. Relationship between the compressive strength (fcu) and the positive ion concentration (C) or pH Value (V). Figure 7 gives the corrosive products for the cemented soil in the H2SO4 solution with a 18.0g/L concentration. As shown, the peaks for the C-S-H (3CaO-2SiO2-3H2O) phase diffraction are small and there is virtually no peak for either Ca(OH)2 (C-H) or 3CaO-Al2O3-3H2O (C-A-H) phase diffraction. For the H2SO4 solution with a high concentration of 18.0g/L, the cementing agent virtually does not function as intended. Based on the above results, the H2SO4 solution is acting in the acid state. These chemical reactions resolve the main cementing agents contained in the cemented soil, i.e., Ca(OH)2, 3CaO-2SiO2-3H2O, and 3CaO-Al2O3-3H2O, and cause the surface peeling and strength reduction. Figure 7. Chemical products in the H2SO4 solution (HS4) with a 18.0g/L concentration and in pure water (W). l:Si02 2:АЬС>з 3:CaSÓ4 4:CaC03 SCaSiOjj 6JSTa2Si03 7:C-S-H S:C-A-S-H JlV Jit wiw Jm 11 ^тл^шятм» i Ž— H S4 »i i.'.1, .tuf i ut im и' Im L ■ Jiii Л ■ i. и - ■ U i Aj w 20 30 40 _ 50 60 Diffraction angle 26/ ° 70 Figure 8. Chemical products in MgSO4 solution (MS4) with 18.0g/L concentration and in pure water (W). Figure 7 also shows the apparently intense peaks for the C-A-S-H and CaSO4 phase diffractions, where C-A-S-H means 3CaO-Al2O3-3CaSO4-32H2O and 3CaO-Al2O3-CaSO4-18H2O. Besides the resolving action, the dissociative SO42- in the solution may take the following chemical reactions. 3CaSO4 + 4CaO-Al2O3-19H2O + 14H2O = (5 3CaO-Al2O3-3CaSO4-32H2O + Ca(OH)2 ( ) 2CaSO4 + 3CaO-Al2O3-CaSO4-18H2O + 14H2O = 3CaO-Al2O3-3CaSO4-32H2O (6) CaCO3 + Ca(OH)2 + SiO2 + CaSO4-2H2O + 12H2O = CaCO3-CaSO4-CaSiO3-15H2O (7) The effects of these reactions on the compressive strength are less significant for the cemented soil in the H2SO4 solution with a lower concentration as the ion H+ plays the main role. Those effects become more significant when the concentration is higher. The crystallizing resultants of 3CaO-Al2O3-3CaSO4-32H2O (C-A-S-H) and CaCO3-CaSO4-CaSiO3-15H2O (C-S-C-H) are apparently larger than the reactant. As such, their inflating force is greater than the sticking force in the cemented soil, and thus causes the cracking in the cemented soil block and the strength reduction. In conclusion, the H2SO4 solution results in a crystallizing corrosion in addition to the resolving corrosion. (2) Corrosion mechanism for the cemented soil in the MgSÜ4 solution Figure 8 gives the corrosive products for the cemented soil in the MgSO4 solution with a 18.0g/L concentration. From the figure, the peaks for the products C-A-S-H, M-A-H (MgO-Al2O3-H2O), M-S-H (MgO-SiO2-H2O), and CaSO4 appear more intense. The main chemical reactions for the cemented soil in the MgSO4 solution are expressed using Equations (8)-(13). 3CaO-2SiO2-3H2O (C-S-H) + 3MgSO4 + 6H2O= 3[CaSO4-2H2O] + 3Mg(OH)2 + 2SiO2 (8) 3CaO-Al2O3-3H2O(C-A-H) + 3MgSO4 + 6H2O= 3[CaSO4-2H2O] + 3Mg(OH)2 + Al2O3 MgCl2 + Ca(OH)2 + 6H2O = CaCl2-6H2O + Mg(OH)2 3Mg(OH)2 + MgCl2 + 8H2O = 2Mg2(OH)3Cl-4H2O 2Mg(OH)2 + 3CaO-2SiO2-3H2O(C-S-H) = 2[MgO-SiO2-H2O] (M-S-H) + 3Ca(OH)2 Mg(OH)2 + 3CaO-Al2O3-3H2O(C-A-H)= MgO-Al2O3-H2O (M-A-H) + 3Ca(OH)2 (9) (10) (11) (12) (13) Equations (8) and (9) indicate a chemical reaction will take place to dissolve the main cementing agent (i.e., 3CaO-2SiO2-3H2O and 3CaO-2Al2O3-3H2O), causing surface peeling for the cemented soil and strength reduction. It has been well documented that the volumes for the resultants CaCl2-6H2O and Mg2(OH)3Cl-4H2O [i.e., Equations (10) and (11)] are seven times those of Ca(OH)2. So, the new products, CaCl2-6H2O and Mg2(OH)3Cl-4H2O, can fill in the voids of the cemented soil. This is beneficial for the unconfined compressive strength of cemented soils when the mass of MgCl2 is suitable and the new reactants (CaCl2-6H2O and Mg2(OH)3Cl-4H2O) are small in each MgSO4 solution. At a 1.5g/L concentration, the volume of CaSO4-2H2O resultants [from Equations (8) and (9)] is two times that of the Ca(OH)2. CaSO4-2H2O is desirable for the compressive strength of cemented soil as the mass of CaSO4-2H2O is deemed suitable. This explains why the strength of the cemented soil cured in MgSO4 solutions with a 1.5g/L concentration is greater. However, the product of Mg(OH)2 may react with 3CaO-2SiO2-3H2O (C-S-H) and 3CaO-Al2O3-3H2O (C-A-H) if there is enough Mg(OH)2 in the MgSO4 solution. The new reactants such as MgO-SiO2-H2O (M-S-H) and MgO-Al2O3-H2O (M-A-H) have poor coagulation, resulting in an unsteady structure and a reduced soil strength. Consequently, the corrosion of the cemented soil in the MgSO4 solution is a resolving corrosion. Besides the resolving action, the dissociative product CaSO4-2H2O [produced by Equations (8) and (9)] may take the chemical reactions as expressed by Equations (5)-(7). The products 3CaO-Al2O3-3CaSO4-32H2O and 3CaO-Al2O3-CaSO4-18H2O (C-A-S-H) and CaSO4-2H2O may also have crystallizing corrosion. So, the corrosion for the cemented soil cured in the MgSO4 solution is a combined resolving and crystallizing corrosion. (3) Corrosion mechanism for the cemented soil in the Na2SO4 solution Figure 9 gives the corrosive products for the cemented soil in the Na2SO4 solution with a 18.0g/L concentration. From the figure, the peaks for products C-A-S-H and CaSO4 are more apparent. In addition to Equations (5)-(7), other chemical reactions for the cemented soil in the Na2S04 solution include Equations (14)-(17). Ut uL.....11 L wf ЦД) l:Si02 2:AL,Ć>3 3:Ca(OH>2 4:CaS04 5:CaC03 6:CaSi03 7:Na2SiÖ, 8C-S-H 9:C-A-S-H • NS4 «J. w 20 30 _ 40 50 Diffraction angle 20/ ' 60 70 Figure 9. Chemical products in Na2SO4 solution (NS4) with 18.0g/L concentration and in pure Water (W). Na2SO4 + 10H2O = Na2SO4-10H2Q Ca(OH)2 + Na2SO4-10H2O = CaSO4-2H2O + 2NaOH + 8H2O 3CaO-2SiO2-3H2O (C-S-H) + 4NaOH = 3Ca(OH)2 + 2Na2SiO3 + 2H2O (14) (15) (16) 3CaO-Al2O3-6H2O(C-A-H)+2NaOH = 3Ca(OH)2 + Na2O-Al2O3 + 4H2O For the Na2SO4 solution with a 1.5g/L concentration, with suitable mass of Na2SO4, the new reactants (i.e., C-A-S-H and CaSO4-2H2O) can fill the voids of the cemented soil and thus increase the unconfined compressive strength. This also explains why the strength of the cemented soil cured in the Na2SO4 solution with a lower concentration at 1.5g/L is larger. For higher concentrations of the Na2SO4 solution, due to the excessive mass of Na2SO4, the Na2SO4 will be combined with H2O to produce the new crystallizing product of Na2SO4-10H2O [Equation (14)]. Na2SO4-10H2O tends to inflate with the cemented soil with C-A-S-H and CaSO4-2H2O, the so-called "crystallizing inflation", and reduces the strength. One the other hand, Na2SO4 combined with Ca(OH)2 produces a new product, NaOH [Equatiion (15)]. NaOH may react with 3CaO-2SiO2-3H2O and 3CaO-Al2O3-6H2O, when the amount of NaOH is enough in the cemented soil [Equations (16) and (17)]. Therefore, the sticking materials in the cemented soil will be decomposed and the new poorly coagulating reactants such as Na2SiO3 and Na2O-Al2O3 will be formed, which dissolve the cemented soil. Hence, the corrosion of the cemented soil cured in the Na2SO4 solution is a combined dissolving and crystallizing corrosion. In summary, the corrosion mechanism for the cemented soil in a H2SO4 or MgSO4 solution is a composite type involving resolving and crystallizing. In addition, it is a composite type consisting of dissolving and crystallizing for the materials in a Na2SO4 solution. 4 CONCLUSIONS_ Based on the results, the following conclusions can be drawn. (1) The sulfate solution changes the cement-soil block including peeling, size reduction, and cracking. The effect of the solution on the block increases with the increase of the concentration of the solution and the curing time. (2) The unconfined compressive strength of the cemen-ted-soil block decreases with the increase of sulfate solution concentration and the curing time. The derived a coefficients can be used to predict the modified compressive strength of cemented soil in various concentrations of corrosive solution. (3) The degree of the corrosion effect from the various solutions is in the descending order: Na2SO4 > MgSO4 > H2SO4. When the SO42- ion exists in a solution, the corrosion degree for the positive ions follows this descending order Na+ > Mg2+ > H+. (4) In terms of the corrosion mechanism for the cemented soil, the corrosion type is found to be a combined resolving and crystallizing for the H2SO4 and MgSO4 solutions and a combined dissolving and crystallizing for the Na2SO4 solution. Acknowledgement This study was financially supported by the National Natural Science Foundation of China (Grant No. 51078253 & 51208333) and Natural Science Foundation of Shanxi Province (2014011036-1, 2014131019, TYUT2014YQ017, 0IT2015). The opinions expressed in this paper are solely of the authors, however. REFERENCES [1] Venkatarama-Reddy, B.V. and Jagadish, K.S. 1995. Influence of soil composition on the strength and durability of soil-cement blocks. Ind. Con. J. 69, 2, 517-524. [2] Walker, P.J. 1995. Strength, durability and shrinkage characteristics of cement stabilized soil blocks. Cement. and Con. Comp.17, 4, 301-310. [3] Pei, X.J.,Yang, G.C. 2000. Research on preventing cement erosion in marine soil. J. Changchun Institute of Tech. 9, 1, 12-14. [4] Shihata, S.A., Baghdadi, Z.A. 2001. Long-term strength and durability of soil cement. J. of Mater. in Civil Eng. 13, 3, 161-165. [5] Jiao, Z.B., Liu, H.L. 2005. Experimental study on cement soil strength in mucky-acid soil. Chin. J. Rock and Soil Mech. 26, 5, 57-60. [6] S. Kolias, V., Kasselouri-Rigopoulou, Karahalios, A. 2005. Stabilization of clayey soils with high calcium fly ash and cement. Cem. Concr. Compos. 27, 2, 301-313. [7] Ning, B.K., Chen, S.L., Liu, B. 2005. Fracturing behaviors of cemented soil under environmental erosion, Chin. J. Rock and Soil Mech. 24, 10, 1778-1782. [8] Ning, B.K., Jin, Chen S.J. 2006. Influence of erosive ions on mechanical properties of cemented soil. J. Shenyang Univ. Technol. 28, 2, 178-181. [9] Dong, X.Q., Bai, X.H., Zhao, Y.Q., Han, P.J. 2007. Study on electrical resistivity of soil-cement polluted by H2SO4. Chin. J. Rock and Soil Mech. 28, 8, 1453-1548. [10] Dong, X.Q., Bai, X.H., Lv, Y.K. 2011. The influence of pH value and SO42- concentration on strength of cemented soil. Adv. Mater. Res. 223, 5, 20062010. [11] Iyengar, S., Al-Tabbaa, A. 2008. Application of two novel magnesia-based cements in the stabilization/ solidification of contaminated soils. Geotech. Spec. Pub. 177, 716-723. [12] Xing, H.F., Yang, X.M., Xu, C. 2000. Strength and microstructure of salt-rich soft soil improved by cement. J. Tongji Univ. 36, 12, 1606-1610. [13] Xing, H.F., Yang, X.M., Xu, C., Ye, G.B. 2009. Strength characteristics and mechanisms of salt-rich soil-cement. Eng. Geol. 103, 1, 33-38. [14] Zandieh, A.R., Yasrobi, S.S. 2010. Study of factors affecting the compressive strength of sandy soil stabilized with polymer. Geotech. and Geologic. Eng. 28, 5, 139-145. [15] Heineck, K.S., Lemos, R.G., Lautenschlager, C.E.R., Consoli, N.C. 2010 Behavior of vertical hydraulic barriers composed by sandy soil, bentonite and cement subjected to alkaline contaminants. Geotech. Spec. Pub. 199, 2462-2471. [16] Voglar, G. E., Lestan, D. 2011. Efficiency modeling of solidification / stabilization of multi-metal contaminated industrial soil using cement and additives. J. of Hazard. Mat. 192, 2, 753-762. [17] Liu, Z.P., Quan, Q.S., Liu, J.W., He, J. 2012. Study on durability of cement soil under brine erosion. Adv. Mat. Res. 446-449, 1, 1858-1863. [18] Yang, Y.Y., Wang, G.H, Xie, S.W. 2012. Effect of magnesium sulfate on the unconfined compressive strength of cement-treated soils. J. of Test. and Eva. 40, 7, 1-8. [19] Yang, Y.Y., Wang, G.H, Xie, S.W., Tu, X.M., Huang, X.G. 2013. Effect of mechanical property of cemented soil under the different pH value. App. Clay Sci. 79, No. 7, 19-24. [20] Gu, B.H., Gao, D.Z., Lin, Z.X., Li, S.Z. 2009. Chinese National Standard GB 50021-2009, Code for investigation of geotechnical engineering. Chinese building industry Press, Beijing. [21] Fan, D.E., He, J.Y., Yang, W.J. 2008. Chinese National Standard GB/T 50046-2008, Code for anticorrosion design of industrial constructions. Chinese Project Press, Beijing. [22] Wang, P.M., Xu, Q.W. 2005. Materials research methods. Chin. Science Press, Beijing. NAVODILA AVTORJEM Vsebina članka Članek naj bo napisan v naslednji obliki: - Naslov, ki primerno opisuje vsebino članka in ne presega 80 znakov. - Izvleček, ki naj bo skrajšana oblika članka in naj ne presega 250 besed. Izvleček mora vsebovati osnove, jedro in cilje raziskave, uporabljeno metodologijo dela, povzetek izidov in osnovne sklepe. - Največ 6 ključnih besed, ki bi morale biti napisane takoj po izvlečku. - Uvod, v katerem naj bo pregled novejšega stanja in zadostne informacije za razumevanje ter pregled izidov dela, predstavljenih v članku. - Teorija. - Eksperimentalni del, ki naj vsebuje podatke o postavitvi preiskusa in metode, uporabljene pri pridobitvi izidov. - Izidi, ki naj bodo jasno prikazani, po potrebi v obliki slik in preglednic. - Razprava, v kateri naj bodo prikazane povezave in posplošitve, uporabljene za pridobitev izidov. Prikazana naj bo tudi pomembnost izidov in primerjava s poprej objavljenimi deli. - Sklepi, v katerih naj bo prikazan en ali več sklepov, ki izhajajo iz izidov in razprave. - Vse navedbe v besedilu morajo biti na koncu zbrane v seznamu literature, in obratno. Dodatne zahteve - Vrstice morajo biti zaporedno oštevilčene. - Predložen članek ne sme imeti več kot 18 strani (brez tabel, legend in literature); velikost črk 12, dvojni razmik med vrsticami. V članek je lahko vključenih največ 10 slik. Isti rezultati so lahko prikazani v tabelah ali na slikah, ne pa na oba načina. - Potrebno je priložiti imena, naslove in elektronske naslove štirih potencialnih recenzentov članka. Urednik ima izključno pravico do odločitve, ali bo te predloge upošteval. Enote in okrajšave V besedilu, preglednicah in slikah uporabljajte le standardne označbe in okrajšave SI. Simbole fizikalnih veličin v besedilu pišite poševno (npr. v, T itn.). Simbole enot, ki so sestavljene iz črk, pa pokončno (npr. Pa, m itn.). Vse okrajšave naj bodo, ko se prvič pojavijo, izpisane v celoti. Slike Slike morajo biti zaporedno oštevilčene in označene, v besedilu in podnaslovu, kot sl. 1, sl. 2 itn. Posnete naj bodo v katerem koli od razširjenih formatov, npr. BMP, JPG, GIF. Za pripravo diagramov in risb priporočamo CDR format (CorelDraw), saj so slike v njem vektorske in jih lahko pri končni obdelavi preprosto povečujemo ali pomanjšujemo. Pri označevanju osi v diagramih, kadar je le mogoče, uporabite označbe veličin (npr. v, T itn.). V diagramih z več krivuljami mora biti vsaka krivulja označena. Pomen oznake mora biti razložen v podnapisu slike. Za vse slike po fotografskih posnetkih je treba priložiti izvirne fotografije ali kakovostno narejen posnetek. Preglednice Preglednice morajo biti zaporedno oštevilčene in označene, v besedilu in podnaslovu, kot preglednica 1, preglednica 2 itn. V preglednicah ne uporabljajte izpisanih imen veličin, ampak samo ustrezne simbole. K fizikalnim količinam, npr. t (pisano poševno), pripišite enote (pisano pokončno) v novo vrsto brez oklepajev. Vse opombe naj bodo označene z uporabo dvignjene številke1. Seznam literature Navedba v besedilu Vsaka navedba, na katero se sklicujete v besedilu, mora biti v seznamu literature (in obratno). Neobjavljeni rezultati in osebne komunikacije se ne priporočajo v seznamu literature, navedejo pa se lahko v besedilu, če je nujno potrebno. Oblika navajanja literature V besedilu: Navedite reference zaporedno po številkah v oglatih oklepajih v skladu z besedilom. Dejanski avtorji so lahko navedeni, vendar mora obvezno biti podana referenčna številka. Primer: ».....kot je razvidno [1,2]. Brandl and Blovsky [4], sta pridobila drugačen rezultat...« V seznamu: Literaturni viri so oštevilčeni po vrstnem redu, kakor se pojavijo v članku. Označimo jih s številkami v oglatih oklepajih. Sklicevanje na objave v revijah: [1] Jelušič, P., Žlender, B. 2013. Soil-nail wall stability analysis using ANFIS. Acta Geotechnica Slovenica 10(1), 61-73. Sklicevanje na knjigo: [2] Šuklje, L. 1969. Rheological aspects of soil mechanics. Wiley-Interscience, London Sklicevanje na poglavje v monografiji: [3] Mitchel, J.K. 1992. Characteristics and mechanisms of clay creep and creep rupture, in N. Guven, R.M. Pollastro (eds.), Clay-Water Interface and Its Rheological Implications, CMS Workshop Lectures, Vol. 4, The clay minerals Society, USA, pp. 212-244.. Sklicevanje na objave v zbornikih konferenc: [4] Brandl, H., Blovsky, S. 2005. Slope stabilization with socket walls using the observational method. Proc. Int. conf. on Soil Mechanics and Geotechnical Engineering, Bratislava, pp. 2485-2488. Sklicevanje na spletne objave: [5] Kot najmanj, je potrebno podati celoten URL. Če so poznani drugi podatki (DOI, imena avtorjev, datumi, sklicevanje na izvorno literaturo), se naj prav tako dodajo. INSTRUCTIONS FOR AUTHORS Format of the paper The paper should have the following structure: - A Title, which adequately describes the content of the paper and should not exceed 80 characters; - An Abstract, which should be viewed as a mini version of the paper and should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation and the methodology employed; it should also summarise the results and state the principal conclusions; - Immediately after the abstract, provide a maximum of 6 keywords; - An Introduction, which should provide a review of recent literature and sufficient background information to allow the results of the paper to be understood and evaluated; - A Theoretical section; - An Experimental section, which should provide details of the experimental set-up and the methods used to obtain the results; - A Results section, which should clearly and concisely present the data, using figures and tables where appropriate; - A Discussion section, which should describe the relationships shown and the generalisations made possible by the results and discuss the significance Podatki o avtorjih Članku priložite tudi podatke o avtorjih: imena, nazive, popolne poštne naslove, številke telefona in faksa, naslove elektronske pošte. Navedite kontaktno osebo. Sprejem člankov in avtorske pravIce Uredništvo si pridržuje pravico do odločanja o sprejemu članka za objavo, strokovno oceno mednarodnih recenzentov in morebitnem predlogu za krajšanje ali izpopolnitev ter terminološke in jezikovne korekture. Z objavo preidejo avtorske pravice na revijo ACTA GEOTECHNICA SLOVENICA. Pri morebitnih kasnejših objavah mora biti AGS navedena kot vir. Vsa nadaljnja pojasnila daje: Uredništvo ACTA GEOTECHNICA SLOVENICA Univerza v Mariboru, Fakulteta za gradbeništvo, prometno inženirstvo in arhitekturo Smetanova ulica 17, 2000 Maribor, Slovenija E-pošta: ags@uni-mb.si of the results, making comparisons with previously published work; - Conclusions, which should present one or more conclusions that have been drawn from the results and subsequent discussion; - A list of References, which comprises all the references cited in the text, and vice versa. Additional Requirements for Manuscripts - Use double line-spacing. - Insert continuous line numbering. - The submitted text of Research Papers should cover no more than 18 pages (without Tables, Legends, and References, style: font size 12, double line spacing). The number of illustrations should not exceed 10. Results may be shown in tables or figures, but not in both of them. - Please submit, with the manuscript, the names, addresses and e-mail addresses of four potential referees. Note that the editor retains the sole right to decide whether or not the suggested reviewers are used. Units and abbreviations Only standard SI symbols and abbreviations should be used in the text, tables and figures. Symbols for physical quantities in the text should be written in Italics (e.g. v, T, etc.). Symbols for units that consist of letters should be in plain text (e.g. Pa, m, etc.). All abbreviations should be spelt out in full on first appearance. Figures Figures must be cited in consecutive numerical order in the text and referred to in both the text and the caption as Fig. 1, Fig. 2, etc. Figures may be saved in any common format, e.g. BMP, JPG, GIF. However, the use of CDR format (CorelDraw) is recommended for graphs and line drawings, since vector images can be easily reduced or enlarged during final processing of the paper. When labelling axes, physical quantities (e.g. v, T, etc.) should be used whenever possible. Multi-curve graphs should have individual curves marked with a symbol; the meaning of the symbol should be explained in the figure caption. Good quality black-and-white photographs or scanned images should be supplied for the illustrations. Tables Tables must be cited in consecutive numerical order in the text and referred to in both the text and the caption as Table 1, Table 2, etc. The use of names for quantities in tables should be avoided if possible: corresponding symbols are preferred. In addition to the physical quantity, e.g. t (in Italics), units (normal text), should be added on a new line without brackets. Any footnotes should be indicated by the use of the superscript1. LIST OF references Citation in text Please ensure that every reference cited in the text is also present in the reference list (and vice versa). Any references cited in the abstract must be given in full. Unpublished results and personal communications are not recommended in the reference list, but may be mentioned in the text, if necessary. Reference style Text: Indicate references by number(s) in square brackets consecutively in line with the text. The actual authors can be referred to, but the reference number(s) must always be given: Example: "... as demonstrated [1,2]. Brandl and Blovsky [4] obtained a different result ..." List: Number the references (numbers in square brackets) in the list in the order in which they appear in the text. Reference to a journal publication: [1] Jelušič, P., Žlender, B. 2013. Soil-nail wall stability analysis using ANFIS. Acta Geotechnica Slovenica 10(1), 61-73. Reference to a book: [2] Šuklje, L. 1969. Rheological aspects of soil mechanics. Wiley-Interscience, London Reference to a chapter in an edited book: [3] Mitchel, J.K. 1992. Characteristics and mechanisms of clay creep and creep rupture, in N. Guven, R.M. Pollastro (eds.), Clay-Water Interface and Its Rheological Implications, CMS Workshop Lectures, Vol. 4, The clay minerals Society, USA, pp. 212-244. Conference proceedings: [4] Brandl, H., Blovsky, S. 2005. Slope stabilization with socket walls using the observational method. Proc. Int. conf. on Soil Mechanics and Geotechni-cal Engineering, Bratislava, pp. 2485-2488. Web references: [5] As a minimum, the full URL should be given and the date when the reference was last accessed. Any further information, if known (DOI, author names, dates, reference to a source publication, etc.), should also be given. Author information The following information about the authors should be enclosed with the paper: names, complete postal addresses, telephone and fax numbers and E-mail addresses. Indicate the name of the corresponding author. Acceptance of papers and copyright The Editorial Committee of the Slovenian Geotechnical Review reserves the right to decide whether a paper is acceptable for publication, to obtain peer reviews for the submitted papers, and if necessary, to require changes in the content, length or language. On publication, copyright for the paper shall pass to the ACTA GEOTECHNICA SLOVENICA. The AGS must be stated as a source in all later publication. For further information contact: Editorial Board ACTA GEOTECHNICA SLOVENICA University of Maribor, Faculty of Civil Engineering, Transportation Engineering and Architecture Smetanova ulica 17, 2000 Maribor, Slovenia E-mail: ags@uni-mb.si NAMEN REVIJE AIMS AND SCOPE Namen revije ACTA GEOTECHNICA SLOVENICA je objavljanje kakovostnih teoretičnih člankov z novih pomembnih področij geomehanike in geotehnike, ki bodo dolgoročno vplivali na temeljne in praktične vidike teh področij. ACTA GEOTECHNICA SLOVENICA objavlja članke s področij: mehanika zemljin in kamnin, inženirska geologija, okoljska geotehnika, geosintetika, geotehnične konstrukcije, numerične in analitične metode, računalniško modeliranje, optimizacija geotehničnih konstrukcij, terenske in laboratorijske preiskave. ACTA GEOTECHNICA SLOVENICA aims to play an important role in publishing high-quality, theoretical papers from important and emerging areas that will have a lasting impact on fundamental and practical aspects of geomechanics and geotechnical engineering. ACTA GEOTECHNICA SLOVENICA publishes papers from the following areas: soil and rock mechanics, engineering geology, environmental geotechnics, geosynthetic, geotechnical structures, numerical and analytical methods, computer modelling, optimization of geotechnical structures, field and laboratory testing. Revija redno izhaja dvakrat letno. The journal is published twice a year. AVTORSKE PRAVICE Ko uredništvo prejme članek v objavo, prosi avtorja(je), da prenese(jo) avtorske pravice za članek na izdajatelja, da bi zagotovili kar se da obsežno razširjanje informacij. Naša revija in posamezni prispevki so zaščiteni z avtorskimi pravicami izdajatelja in zanje veljajo naslednji pogoji: Fotokopiranje V skladu z našimi zakoni o zaščiti avtorskih pravic je dovoljeno narediti eno kopijo posameznega članka za osebno uporabo. Za naslednje fotokopije, vključno z večkratnim fotokopiranjem, sistematičnim fotokopiranjem, kopiranjem za reklamne ali predstavitvene namene, nadaljnjo prodajo in vsemi oblikami nedobičk-onosne uporabe je treba pridobiti dovoljenje izdajatelja in plačati določen znesek. Naročniki revije smejo kopirati kazalo z vsebino revije ali pripraviti seznam člankov z izvlečki za rabo v svojih ustanovah. COPYRIGHT Upon acceptance of an article by the Editorial Board, the author(s) will be asked to transfer copyright for the article to the publisher. This transfer will ensure the widest possible dissemination of information. This review and the individual contributions contained in it are protected by publisher's copyright, and the following terms and conditions apply to their use: Photocopying Single photocopies of single articles may be made for personal use, as allowed by national copyright laws. Permission of the publisher and payment of a fee are required for all other photocopying, including multiple or systematic copying, copying for advertising or promotional purposes, resale, and all forms of document delivery. Subscribers may reproduce tables of contents or prepare lists of papers, including abstracts for internal circulation, within their institutions. Elektronsko shranjevanje Za elektronsko shranjevanje vsakršnega gradiva iz revije, vključno z vsemi članki ali deli članka, je potrebno dovoljenje izdajatelja. Electronic Storage Permission of the publisher is required to store electronically any material contained in this review, including any paper or part of the paper. ODGOVORNOST Revija ne prevzame nobene odgovornosti za poškodbe in/ali škodo na osebah in na lastnini na podlagi odgovornosti za izdelke, zaradi malomarnosti ali drugače, ali zaradi uporabe kakršnekoli metode, izdelka, navodil ali zamisli, ki so opisani v njej. RESPONSIBILITY No responsibility is assumed by the publisher for any injury and/or damage to persons or property as a matter of product liability, negligence or otherwise, or from any use or operation of any methods, products, instructions or ideas contained in the material herein. .л О ч<2> .МЛ/ .-A/ v ^ VT \ 4- è .S* (В cS "V«»*0 SC / iT «V ^ y <4« 9771854017001