Numerical Simulation of Ground Water Flow in Dual Porous Media of the Karun 4 Dam (Iran) Foundation and Abutments Dissertation for the acquisition of the academic degree Doktor der Ingenieurswissenschaften (Dr.-Ing.) Submitted to the Faculty of Civil and Environmental Engineering of the University of Kassel by Seyed Mohammad Hosseiny Sohi First supervisor: Prof. Dr. rer. nat. Manfred Koch Second supervisor: Assoc. Prof. Dr. Javad Ashjari Defense Date: 30.10.2018 Kassel, Germany March 2019 I Erklärung Hiermit versichere ich, dass ich die vorliegende Dissertation selbständig, ohne unerlaubte Hilfe Dritter angefertigt und andere als die in der Dissertation angegebenen Hilfsmittel nicht benutzt habe. Alle Stellen, die wörtlich oder sinngemäß aus veröffentlichten oder unveröffentlichten Schriften entnommen sind, habe ich als solche kenntlich gemacht. Dritte waren an der inhaltlichen Erstellung der Dissertation nicht beteiligt; insbesondere habe ich nicht die Hilfe eines kommerziellen Promotionsberaters in Anspruch genommen. Kein Teil dieser Arbeit ist in einem anderen Promotions- oder Habilitationsverfahrendurch mich verwendet worden. Seyed Mohammad Hosseiny Sohi Kassel, October 2018 II Acknowledgements I would like to extend thanks to the many people, from many countries, Germany, Iran, Turkey, Colombia, Pakistan and USA, who so generously contributed to the work presented in this thesis. Firstly, I would like to express my sincere gratitude to my advisor Prof. Manfred Koch for the continuous support of my Ph.D. study and related research, for his patience, motivation, and immense knowledge. His guidance helped me in all the time of research and writing of this thesis. I am also grateful to my second advisor, Dr. Javad Ashjari, from the University of Kharazmi, Tehran, Iran. I am extremely thankful and indebted to him for sharing expertise, and sincere and valuable guidance and encouragement extended to me. I gratefully acknowledge the funding sources that made my Ph.D. work possible. I was funded by the Rosa Luxembourg foundation for my first 3 years. Many thanks to Prof. Aram Ziai the liaison lecturer of the foundation at the university of Kassel, who interviewed me by application for the scholarship. I’m extremely grateful to Dr. Sandra Thieme, the responsible of STEM disciplines of the foundation, who helped continuously during the promotion. Very strong thanks to Mr. Etaati from Iran Water and Power resources development Co (IWPC) provided me the access to the data and reports. Special Thanks to Dr. Richard Winston from U.S. Geological Survey, the developer of the ModelMuse. I take this opportunity to express gratitude to all the department faculty members for their help and support. The secretary Mrs. Astrid Inden was always ready to help and support the official procedures. Thanks to Mrs. Monika Pruditsch my office was always clean and comfortable. Dr. Majid Taei Semiromi was always ready for constructive scientific discussions and my endless questions. Tim Wolters with his great open-minded behavior and his support left a very nice image of a kind colleague in my mind. Maria with her power of life was the reason of III resistance. Yakup Akcuru and his nice smiles was making my days at the mornings. Shiva Mirdar with her cookies brought the sweetness to our department. I appreciate Mr. Bijan Ferdowsi who was in Kassel the meaning of having a reliable friend, in all circumstances. A special thanks to my family. Words cannot express how grateful I am to my mother, Fahim (although she is gone), and father, Mostafa, my mother-in law, Najmeh, father-in-law, Hamid, my sister Zahra and my brothers Barman, Masoud, Mohammad Hossein for all the sacrifices that you’ve made on my behalf. And finally, this work will be dedicated to my best friend and beloved wife, Maryam Soleimani Rad who taught me that with love, a lot of impossibilities become possible. Seyed Mohammad Hosseiny Sohi IV Abstract When it comes to the construction of dams and reservoirs in karsts, which are known for their high perviousness due to dissolution faults and conduits, the prevailing risk is water loss. Water loss through basements and abutments of the dams constructed in karst areas leads to considerable costs and non-productivity of the dams. Although the construction of a grout curtain is the most common method to reduce or prevent the water loss, it can solve the problem only partly and in some cases this operation is even fruitless. In this regard, a numerical model, which simulates the dual porous media of carbonate limestone, could prognose the possible water leakage paths. In the present study, different stages to reach this purpose are described for a dam site in Iran as follows: 1) The Karun 4 concrete dam is with a maximum dam height above foundation of 230m the highest dam of Iran. It has been in operation since July 2011. It was constructed over a karst limestone and water leakage has been, and continues to be, one of the most important problems of this dam. The geological formations around the dam site consist of carbonate layers of the Asmari formation and marlstone and marly limestone of the impervious Pabdeh formation, upstream of the dam site. The object of this study is to determine the dam construction effects on the groundwater levels in the adjacent aquifer. Thus, before the dam construction, the natural groundwater levels at both banks were recorded to be 5-8m higher than the river stage. This indicates that the natural groundwater gradient at that time was from the banks towards the river, i.e. the latter was an effluent stream before the construction of dam. At that time no important springs were identified at the dam site., other than a minor spring with a discharge <5 lit/sec that emanated at the contact interface between the permeable Asmari- and the impermeable Pabdeh formation (perching aquifer conditions) on the left bank, upstream of the dam. A Lugeon test that was carried out at that time indicated that the permeability of the adjacent limestone is high, as it varies from 25 to 55 Lu in the layers above the river level, but it decreases gradually to 3 Lu in the formations below. After dam impounding, some changes in the borehole's water levels were observed. Thus, it has been found, in particular, that the leakage from the reservoir has induced groundwater level rises between 12 to 17 meters. 2) Karun 4 dam is the highest concrete double arc dam in Iran located in the Zagros mountain range in southwest of Iran. It has been under operation since March 2011. The dam foundations and abutments are within the Asmari formation, consisting of fractured carbonate rocks with a high potential of karstification. The water tightness of the dam was one of the main issues during the construction. The grout curtain was constructed to solve the problem in three steps. The object of this study is to evaluate the designed curtain based on the grout taking rate, permeability test, and geological logs. The studied boreholes consist firstly of the exploration boreholes, which represent the structural situation and characteristics of the rock formation, and, secondly of the check holes that were drilled after finalizing the implementation of the grout curtain. The results illustrate the efficiency of the grout curtain, but, nevertheless, some water leakage is still occurring through the curtain, so that further monitoring and treatment is required to ensure the dam safety, namely, to prevent the water leakage. 3) Karun 4 concrete dam with 230 m height is the first dam constructed at the upstream of Karun river in south west of Iran. The dam site area is located on the carbonate limestone of the Asmari formation which has strong karstic characteristics and, subsequently, a high potential for water leakage. To deal with this problem, a seepage controlling system is implemented within the dam that includes a grout curtain as a virtual barrier against water leakage on the dam foundation and abutments, and a drainage curtain at the downstream of the dam body to allow collection of the escaped water from the grout curtain; and the monitoring instruments: Standpipe piezometers and observation wells. The object of present study is to present the dam seepage controlling system; and to investigate the water head variation of the piezometers and observation wells, using obtained measurements after impounding, to evaluate the efficiency of the grout curtain. In addition, the areas which have larger water head variations are discussed, indicating the variations are generally higher in the left- than in the right bank. And finally, the need of a numerical model of the groundwater flow of the whole area to project the dam reservoir behavior for the coming years, is concluded. 4) Water leakage through dam foundation and abutments is one of the prevalent problems when a dam is constructed in a carbonate aquifer with karst potential. The Karun 4 dam in Iran is built in a porous limestone formation that could lead to water leakage. The geological field investigations before its construction indicated the zones with fractures and conduits. To mitigate this problem, a large grout curtain has been drilled and injected with cement grout via different grouting galleries in the dam foundation and abutments. To evaluate the function of the curtain, and find out the possible water leakage paths, a numerical model of dual porous media is set up using Conduit Flow Process (CFP) for MODFLOW-2005. The model simulates both porous matrix laminar- and conduits laminar or turbulent flow. Using the reservoir water levels, drainage observations, observed groundwater head from installed piezometers and observation wells, the model is successfully calibrated and validated. In the transient state, firstly, the matrix medium is calibrated in MODFLOW as a laminar groundwater flow; then, the conduits are set in the CFP model. The conduit layout is determined by examination of the faults from geological maps. The approved faults are assigned as pipes in the model. The pipes calibration parameters are their hydraulic conductivity and diameter. Afterwards the efficiency of the grout curtain is evaluated which identifies various hydraulic conductivities along the grout curtain. The present model enables one to investigate different scenarios of water leakage in the model area which is the subject of a subsequent study. 5) Despite all recent achievements in Karst aquifers’ research, the numerical modelling of groundwater flow processes in such aquifers has remained complicated and problematic, since the precise determination of the conduits network within a particular karst aquifer requires a significant amount of local data, observations by site experts, and sound application findings from laboratories/academic institutes researchers. The objective of the present study is to investigate different scenarios of water leakage through the karstic dam foundation and abutments of Karun 4 dam in Iran. To this avail, classical groundwater flow in the porous matrix medium and laminar or turbulent pipe flow in the discrete conduits are simulated with the recently developed numerical groundwater flow model MODFLOW-CFP. A total of 97 scenarios, divided into three categories, are simulated and the computed amount of downstream water escape compared with that of the calibrated transient reference model scenario (RS): In the first category, for each of the 15 originally recognized karstic faults, represented by calibrated pipes in the reference CFP-model, four different diameters: 20, 50, 200 and 500 % of the corresponding reference diameter of the pipe, are analyzed. In the second category, 10 new hypothetical pipes are embedded at the location of extra faults with three scenarios for each pipe. In the third category several adjusted pair-pipes that represent continuous faults on the up- and downstream of the grout curtain are connected to each other. Results of the first group of scenarios identified the most critical fault, namely, F4-Up, whose pipe characteristics have the largest effect on the downstream water escape. Results of the second scenario category emphasize that some potential faults at the border of the adjacent As1-g and As1-i karstic formation units may lead to an additional water discharge of nearly 1% relative to that of the RS. Finally, the connected-pipes scenarios results demonstrate the effectiveness of the grout curtain in impeding downstream water leakage even if potential karstic faults close to the critical fault F4-Up may open up further. V Kurzfassung Beim Bau von Staudämmen und Reservoirs in Karst, die für ihre hohe Pererverenz aufgrund von Auflösungsfehlern und Leitungen bekannt sind, besteht das vorherrschende Risiko im Wasserverlust. Der Wasserverlust durch das Fundament und die Widerlager der in den Karstgebieten errichteten Dämme führt zu erheblichen Kosten und zur Nichtproduktivität der Staudämme. Obwohl die Konstruktion eines Dichtungsschleiers die gebräuchlichste Methode zur Verringerung oder Verhinderung von Wasserverlust ist, kann sie das Problem nur teilweise lösen, und in einigen Fällen ist dieser Vorgang sogar fruchtlos. In diesem Zusammenhang könnte ein numerisches Modell, das die dualen porösen Medien von Kalkstein simuliert, mögliche Wasseraustritt Wege vorhersagen. In dieser Studie werden die verschiedenen Stufen, um dieses Ziel zu erreichen, für einen Dammstandort im Iran wie folgt beschrieben: 1) Der Karun 4-Staudamm ist mit einer maximalen Staudammhöhe von über 230 Meter der höchste Staudamm Irans. Er ist seit Juli 2011 in Betrieb. Er wurde über einem karstigen Kalkstein gebaut und die Undichtigkeit zählt zu den größten Problemen dieses Staudamms. Die geologischen Formationen um den Damm herum bestehen aus Karbonatschichten der Asmari-Formation und Mergel(kalk)kalk der undurchlässigen Pabdeh-Formation flussaufwärts des Damms. Ziel dieser Studie ist es, die Auswirkungen der Staudämme auf die Grundwasserspiegel im angrenzenden Grundwasserleiter zu ermitteln. So wurden vor dem Staudammbau die natürlichen Grundwasserstände an beiden Ufern 5 bis8 m höher als in der Flussphase gemessen. Dies deutet darauf hin, dass der natürliche Grundwassergradient zu dieser Zeit von den Ufern in Richtung des Flusses gerichtet war, d. h. letzterer war ein effluenter Fluss vor dem Bau des Staudamms. Zu dieser Zeit wurden keine wichtigen Quellen an der Dammstelle identifiziert, mit Ausnahme einer kleinen Quelle mit einer Austritt von weniger als fünf Litern pro Sekunde, die an der Kontaktfläche zwischen der durchlässigen Asmari - und der undurchlässigen Pabdeh - Formation (schwebender Grundwasserleiter) am linken Ufer, stromaufwärts des Dammes, auftrat. Ein Lugeon-Test, der damals durchgeführt wurde, zeigte, dass die Permeabilität des angrenzenden Kalksteins hoch ist, indem er in den Schichten über dem Fluss von 25 bis 55 Lu variiert, aber in den darunter liegenden Formationen allmählich auf 3 Lu abfällt. Nach dem Aufstauen des Staudamms wurden einige Veränderungen in den Wasserständen des Bohrlochs beobachtet. So zeigt sich insbesondere, dass durch die Undichtigkeit aus dem Reservoir die Grundwasserstände um 12 bis 17 Meter ansteigen. 2) Der Karun 4-Damm ist der höchste Beton-Doppelbogen-Damm im Iran, der sich im Zagros-Gebirge im Südwesten des Iran befindet. Es ist seit März 2011 in Betrieb. Die Dammfundamente und Widerlager liegen innerhalb der Asmari-Formation, die aus gebrochenen Karbonatgesteinen mit hohem Verkarstungspotential besteht. Die Undichtigkeit des Staudamms war eines der Hauptprobleme während des Baus. Der Dichtungsschleier wurde konstruiert, um das Problem in drei Schritten zu lösen. Das Ziel dieser Studie ist es, den entworfenen Dichtungsschleier basierend auf der Rate der Fugenaufnahme, dem Durchlässigkeitstest und den geologischen Logs zu bewerten. Die untersuchten Bohrlöcher bestehen zum einen aus den Erkundungsbohrungen, die die strukturelle Situation und Merkmale der Gesteinsformation darstellen, und zum anderen aus den Kontrollbohrungen, die nach Abschluss der Ausführung des Dichtungsschleiers gebohrt wurden. Die Ergebnisse veranschaulichen die Effizienz des Dichtungsschleiers, dennoch tritt ein gewisser Wasseraustritt durch den Gürtel auf, so dass eine weitere Überwachung und Behandlung erforderlich ist, um die Dammsicherheit zu gewährleisten, also den Wasseraustritt zu verhindern. 3) Der Karun 4-Staudamm ist mit 230 Meter Höhe der erste Damm, der am Fluss Karun im Südwesten des Iran gebaut wurde. Das Dammstandort befindet sich auf dem Karbonatkalk der Asmari-Formation, der starke Karsteigenschaften und folglich ein hohes Potential zur Undichtigkeit aufweist. Um dieses Problem zu lösen, wird ein System zur Kontrolle der Versickerung innerhalb des Staudamms implementiert, das einen Dichtungsschleier als virtuelle Barriere gegen Wasseraustritt auf dem Dammfundament und den Widerlagern und einen Sickerleitung stromabwärts des Dammkörpers enthält, um die Sammlung vom austretenden Wasser aus dem Dichtungsschleier zu ermöglichen; und die Überwachungsinstrumente: Standrohr-Piezometer und Beobachtungsbrunnen. Das Ziel der vorliegenden Studie ist es, das Dammversickerungssteuerungssystem zu präsentieren und die Änderungen der Standrohrspiegelhöhen der Piezometer und Beobachtungsbrunnen unter Verwendung von erhaltenen Messungen nach dem Aufstauen zu untersuchen, um die Effizienz des Mörtelvorhangs zu bewerten. Zusätzlich werden die Bereiche, die größere Änderungen der Standrohrspiegelhöhen aufweisen, diskutiert. Im Allgemeinen sind die Variationen am linken höher als am rechten Ufer. Schließlich wird die Notwendigkeit eines Grundwasserströmungsmodells des gesamten Gebiets zur Darstellung des Dammreservoir-Verhaltens für die kommenden Jahre herausgestellt. 4) Wasseraustritt durch das Dammfundament und die Widerlager ist eines der vorherrschenden Probleme, wenn ein Damm in einem Karbonat-Aquifer mit Verkarstungspotential errichtet wird. Der Karun-4-Damm im Iran ist in einer porösen Kalksteinformation gebaut, die zu Wasseraustritt führen könnte. Die geologischen Felduntersuchungen vor dem Bau zeigten die Zonen mit Brüchen und Gängen/Kanälen. Um dieses Problem zu lindern, wurde ein großer Dichtungsschleier gebohrt und mit Zementmörtel über verschiedene Injektionsgalerien im Dammfundament und den Widerlagern injiziert. Um die Funktion des Dichtungsschleiers zu bewerten und mögliche Wasserleckpfade herauszufinden, wird ein numerisches Modell von zwei porösen Medien unter Verwendung des Conduit Flow Process (CFP) für MODFLOW-2005 aufgebaut. Das Modell simuliert sowohl die laminare Strömung in der porösen Matrix als auch die laminare oder turbulente Strömung in den Kanälen. Unter Verwendung des Wasserpegels des Reservoirs, der Entwässerungsbeobachtungen, des beobachteten Grundwasserstandes von installierten Piezometern und Beobachtungsbrunnen wird das Modell erfolgreich kalibriert und validiert. Im transienten Zustand wird zunächst das Matrixmedium in MODFLOW als laminare Grundwasserströmung kalibriert, danach werden die Kanäle im CFP-Modell festgelegt. Das Kanalnetz wird durch Untersuchungen der Verwerfungen aus geologischen Karten bestimmt. Die so bestätigten Verwerfungen werden im Modell als „Pipes“ (Rohre) zugewiesen. Die Parameter für die Kalibrierung der Rohre sind deren hydraulische Leitfähigkeit und Durchmesser. Anschließend wird der Wirkungsgrad des Dichtungsschleiers bewertet, wodurch verschiedene hydraulische Leitfähigkeiten entlang des Dichtungsschleiers identifiziert werden. Das vorliegende Modell ermöglicht es, verschiedene Szenarien von Wasseraustritten im Modellgebiet zu untersuchen, die Gegenstand der nachfolgenden Studie sind. 5) Trotz aller neueren Erkenntnisse in der Erforschung von Karstgrundwasserleitern ist die numerische Modellierung von Grundwasserströmungsprozessen in solchen Grundwasserleitern sehr kompliziert und problematisch geblieben, da die genaue Bestimmung des Karst- Kanalnetzes in einem bestimmten Karstgrundwasserleiter eine erhebliche Menge lokaler Daten und vertiefte Kenntnisse, sowohl von Experten vor Ort, als auch von Forschern in Laboren / akademischen Instituten, erfordert. Das Ziel der vorliegenden Studie ist es, verschiedene Szenarien des Wasseraustritts durch das Karst- Fundament und -Widerlager des Karun 4-Staudamms im Iran zu untersuchen. Zu diesem Zweck wird die klassische Grundwasserströmung im porösen Matrixmedium und die laminare oder turbulente Rohrströmung in den diskreten Kanälen mit dem kürzlich entwickelten numerischen Grundwasserströmungsmodell MODFLOW-CFP simuliert. Insgesamt 97 Szenarien, unterteilt in drei Kategorien, werden simuliert und die berechnete Menge des nachgelagerten Wasserabflusses mit der des kalibrierten transienten Referenzmodells (RS) verglichen: In der ersten Kategorie werden für jeden der 15 ursprünglich erkannten Karstfugen (Kanälen), die in dem Referenz-CFP-Modell durch kalibrierte Rohre („Pipe“) dargestellt sind, vier verschiedene Durchmesser analysiert: 20, 50, 200 und 500% des entsprechenden Referenzdurchmessers des Rohrs. In der zweiten Kategorie sind 10 neue hypothetische Rohre an Stellen zusätzlicher Karstkanälen mit drei Szenarien für jedes Rohr eingebettet worden. In der dritten Kategorie sind mehrere angepasste Paarrohre miteinander verbunden, die fortlaufende Karstspalten vor und hinter dem Dichtungsschleier darstellen sollen. Die Ergebnisse der ersten Gruppe von Szenarien identifizierten die kritischste Karstfuge, nämlich F4-Up, deren Rohrleitungseigenschaften den größten Einfluss auf den nachgelagerten Wasseraustritt haben. Die Ergebnisse der zweiten Szenario-Kategorie betonen, dass einige potenzielle Karstspalten an der Grenze der benachbarten Karstformationseinheiten As1-g und As1-i zu einem zusätzlichen Wasserabfluss von fast 1% relativ zu dem des RS führen können. Schließlich demonstrieren die Ergebnisse der dritten Szenarien-Gruppe mit verbundenen Rohren, dass der Dichtungsschleier die Wasserleckage nachgeschaltet verhindert, selbst wenn potenzielle Karstspalten in der Nähe der kritischen Spalte F4-Up sich weiter öffnen könnten. VI TABLE OF CONTENTS Acknowledgements .............................................................................................................................. II Abstract ................................................................................................................................................ IV Kurzfassung .......................................................................................................................................... V List of Figures ................................................................................................................................... X List of Tables ................................................................................................................................. XII General Introduction ........................................................................................................................ 1 I. Background ................................................................................................................................. 1 II. Case study: Karun 4 dam / Iran ................................................................................................... 2 III. Research objectives ................................................................................................................. 2 IV. Structure of dissertation .......................................................................................................... 3 Chapter 1- Construction Effects of the Karun 4 Dam, Iran, on the Groundwater in the adjacent Karstic Aquifer .................................................................................................................. 4 1.1 Introduction ............................................................................................................................... 5 1.2 Study area .................................................................................................................................. 6 1.2.1 Site specification and geology ............................................................................................ 6 1.2.2 Hydrogeology ..................................................................................................................... 9 1.3 Study method ............................................................................................................................. 9 1.4 Results ...................................................................................................................................... 10 1.4.1 Before water impound ....................................................................................................... 10 1.4.2 During and after water impound ....................................................................................... 11 1.4.3 Hydrochemistry of the groundwater ................................................................................. 12 1.5 Conclusions and outlook ......................................................................................................... 13 References ............................................................................................................................................ 15 Chapter 2- Evaluating Permeability and Groutability at the Karun 4 Dam, Iran, using Lugeon Values and Grout Take ..................................................................................................... 17 2.1 Introduction ............................................................................................................................. 18 2.2 Study area and dam specifications ........................................................................................ 19 2.2.1 Site specification and geology .......................................................................................... 19 2.2.2 Dam specifications and grout curtain ................................................................................ 20 2.3 Methodology and data ............................................................................................................ 22 2.3.1 Water pressure test (WPT) or Lugeon test ........................................................................ 22 2.3.2 Take ................................................................................................................................... 22 2.3.3 Ewert method .................................................................................................................... 22 2.3.4 Data sources ...................................................................................................................... 22 2.4 Results ...................................................................................................................................... 23 2.4.1 Lugeon permeability ......................................................................................................... 23 VII 2.4.2 Groutability ....................................................................................................................... 25 2.4.3 Correlation of Lu- permeability and Take- groutability based on Ewert’s method .......... 26 2.5 Conclusions .............................................................................................................................. 29 References ............................................................................................................................................ 30 Chapter 3- Monitoring of Seepage Controlling System of Karun 4 Dam in Iran ..................... 33 3.1 Introduction ............................................................................................................................. 34 3.2 Seepage controlling system of the Karun 4 dam .................................................................. 35 3.2.1 Site specification and geology .......................................................................................... 35 3.2.2 Grout curtain ..................................................................................................................... 37 3.2.3 Drainage curtain ................................................................................................................ 38 3.2.4 Monitoring instruments ..................................................................................................... 39 3.2.4.1 Standpipe piezometers .............................................................................................. 39 3.2.4.2 Observation wells ..................................................................................................... 40 3.3 Data .......................................................................................................................................... 40 3.4 Results ...................................................................................................................................... 40 3.4.1 Standpipe piezometers ...................................................................................................... 40 3.4.2 Observation wells .............................................................................................................. 44 3.5 Conclusions .............................................................................................................................. 46 References ............................................................................................................................................ 47 Chapter 4 - Simulation of Groundwater Flow in the Dual Porous Media of the Karun 4 Dam Foundation and Abutments using Conduit Flow Process (CFP) for MODFLOW-2005.......... 49 4.1 Introduction ............................................................................................................................. 50 4.2 Study area ................................................................................................................................ 53 4.2.1 Site specifications and geological setting ......................................................................... 53 4.2.2 Seepage controlling system ............................................................................................... 54 4.3 Numerical model ..................................................................................................................... 56 4.3.1 Conduit Flow Process (CFP) ............................................................................................. 56 4.3.2 Governing equation for porous medium ........................................................................... 56 4.3.3 Governing equations for pipe flow ................................................................................... 58 4.3.4 Graphical User Interface ModelMuse ............................................................................... 60 4.3.5 Set-up and designing of the numerical model ................................................................... 60 4.3.5.1 Conceptual model ..................................................................................................... 60 4.3.5.2 Model grid ................................................................................................................ 60 4.3.5.3 Boundary and initial conditions ................................................................................ 63 4.3.5.4 Hydraulic conductivity ............................................................................................. 63 4.3.5.5 Drainage conductance ............................................................................................... 64 4.3.5.6 Head observation package (HOB) ............................................................................ 66 4.3.5.7 Storage ...................................................................................................................... 67 VIII 4.4 Model calibration .................................................................................................................... 68 4.4.1 Calibration parameters of porous medium in steady state ................................................ 68 4.4.2 Calibration parameters of porous medium in transient state ............................................. 68 4.4.3 CFP calibration (transient state) ........................................................................................ 69 4.5 Results ...................................................................................................................................... 72 4.5.1 Steady state calibration ..................................................................................................... 72 4.5.1.1 Hydraulic conductivity (K) ....................................................................................... 72 4.5.1.2 Drain conductance (C) .............................................................................................. 73 4.5.1.3 Groundwater head ..................................................................................................... 73 4.5.2 Transient state calibration ................................................................................................. 75 4.5.2.1 Storage ...................................................................................................................... 75 4.5.2.2 Grout-curtain-efficiency evaluation .......................................................................... 76 4.5.2.3 Groundwater heads ................................................................................................... 77 4.5.3 Impact of faults on the groundwater leakage through the matrix ..................................... 82 4.6 Conclusions .............................................................................................................................. 87 References ............................................................................................................................................ 88 Chapter 5 - Investigation of different Scenarios of Water Leakage through Foundation and Abutments of Karun 4 Dam ........................................................................................................... 92 5.1 Introduction ............................................................................................................................. 93 5.2 Study area ................................................................................................................................ 96 5.2.1 Site specifications and geological setting ......................................................................... 96 5.2.2 Seepage controlling system ............................................................................................... 98 5.3 Mathematical and numerical methodology .......................................................................... 99 5.3.1 Theoretical basics of the MODFLOW-CPF dual porosity model ..................................... 99 5.3.2 Set-up of the CPF- model ................................................................................................ 102 5.3.2.1 General layout and boundary conditions ................................................................ 102 5.3.2.2 Incorporation of karst conduits and final calibration of the CPF-model ................ 104 5.4 Simulation scenarios of pipe adjustments for the representation of faults ...................... 108 5.4.1 Overview of scenario specifications and definition of the reference scenario ................ 108 5.4.2 Scenarios with changing pipe diameters ......................................................................... 109 5.4.2.1 Pipe discharge results ............................................................................................. 109 5.4.2.2 Drainage results ...................................................................................................... 112 5.4.2.3 Discharges at downstream boundaries (Constant Head OUT) ............................... 112 5.4.2.4 Storage .................................................................................................................... 112 5.4.2.5 Discussion of pipe diameter variation scenarios ..................................................... 112 5.4.3 Scenarios with hypothetical pipe locations to represent possible critical faults ............. 114 5.4.3.1 Model- and pipe specifications ............................................................................... 114 5.4.3.2 Results and discussions........................................................................................... 115 IX 5.4.4 Scenarios with connected pipes (continuous faults) through the grout curtain ............... 119 5.4.4.1 Rational and approach ............................................................................................ 119 5.4.4.2 Results .................................................................................................................... 120 5.5. Conclusions ............................................................................................................................ 122 References .......................................................................................................................................... 124 X List of Figures Figure 1.1 Location of Karun 4 dam in Iran .......................................................................................... 7 Figure 1.2 Karun 4 dam geological plan ................................................................................................ 8 Figure 1.3 Karun 4 dam geological cross-section .................................................................................. 8 Figure 1.4 Change of water levels after water impounding of Karun 4 dam (2010) ........................... 12 Figure 1.5 Variations of hydrochemical parameters between 2010 and 2014 ..................................... 14 Figure 2.1 Location and details of the Karun 4 dam in Iran ................................................................ 20 Figure 2.2 Selected Exploratory- , line 1 and line 2 check holes ......................................................... 21 Figure 2.3 Ewert graph of water pressure test values (Lugeon) over Take .......................................... 23 Figure 2.4 Frequency distribution of Take in the different classes ...................................................... 25 Figure 2.5 Frequency distribution of Take in the different classes ...................................................... 26 Figure 2.6 Lugeon versus Take with the four Ewert’s groupings ........................................................ 27 Figure 2.7 Frequency distribution of Lu-Take pairs in the different groups ........................................ 28 Figure 3.1 Location of Karun 4 dam in Iran with observation wells’ layout ....................................... 36 Figure 3.2 Grout curtain expanded section view and the piezometers’ layout .................................... 37 Figure 3.3 Cross-sectional view of the grout curtain system ............................................................... 39 Figure 3.4 Water heads measured by the downstream piezometers..................................................... 42 Figure 3.5 Water head variations of the four piezometers of each gallery in left and right bank ........ 43 Figure 3.6 Groundwater level measurements of the observation wells ............................................... 45 Figure 4.1 Plan view of the Karun 4 dam and grout curtain in Iran ..................................................... 53 Figure 4.2 Grout curtain expanded cross section ................................................................................. 55 Figure 4.3 Darcy’s law for a prism of porous material (from McDonald and Harbaugh, 1988) ......... 57 Figure 4.4 Conceptual model of Karun 4 dam groundwater flow and boundary conditions layout .... 61 Figure 4.5 Plan view of the 4 layers of the Karun 4 MODFLOW ....................................................... 62 Figure 4.6 Reservoir water level hydrograph from time of impounding till end of modelling time .... 63 Figure 4.7 Geological map for layer 2 including the Asmari 1 subunits ............................................. 64 Figure 4.8 Buried drain pipe (from McDonald and Harbaugh, 1988) ................................................. 65 Figure 4.9 The calibrated zoning sections of the drainage curtain ....................................................... 66 Figure 4.10 Layout of piezometers and observation wells ................................................................... 67 Figure 4.11 Karun4 dam area plan view with faults sketched in red ................................................... 70 Figure 4.12 Layout of the pipes at a) layer 2, b) layer 3, and c) layer 4, and d) faults plan view ........ 72 Figure 4.13 Calibrated K-values of in a) layer 1, b) layer 2, c) layer 3, and, d) layer 4. ..................... 74 Figure 4.14 Steady-state calibrated- versus observed heads (left panel) and residuals (right panel) . 75 Figure 4.15 Steady-state calibrated ground water heads in layer 3 ...................................................... 75 Figure 4.16 Calibrated Ss-values of in the (confined) a) layer 2, b) layer 3, c) layer 4. ...................... 76 Figure 4.17 Observed- and simulated groundwater head time series ................................................... 78 XI Figure 4.18 Residual of observed- and simulated groundwater head ................................................. 79 Figure 4.19 Simulated- versus observed head values (left) and residuals of each observation point (right) – Transient state (calibration phase) .......................................................................................... 80 Figure 4.20 Similar to Figure 4.19, but for the validation phase ......................................................... 80 Figure 4.21 Similar to Figure 4.19 (left), but for well OW1 and OW4 for the calibration phase ....... 80 Figure 4.22 Simulated groundwater head contours of the last stress period on March 13th 2017 in : a) layer 2, b) layer 3 and c) layer 4. .......................................................................................................... 81 Figure 4.23 Water heads in the nodes of the various pipes and the porous matrix medium, and the exchange flow between them, on March 13th 2017, the last stress period of the model ...................... 83 Figure 4.23 (continued) ....................................................................................................................... 84 Figure 4.23 (continued) ........................................................................................................................ 85 Figure 5.1 Plan view of the Karun 4 dam in Iran ................................................................................. 97 Figure 5.2 Grout curtain expanded cross section ................................................................................. 98 Figure 5.3 Darcy’s law for a prism of porous material (from McDonald and Harbaugh, 1988) ....... 100 Figure 5.4 Boundary conditions layout of the numerical model ........................................................ 102 Figure 5.5 Geological map at layer 2 including Asmari 1 subunits ................................................... 103 Figure 5.6 Layout of piezometers and observation wells ................................................................... 104 Figure 5.7 Layout of the pipes at a) layer 2, b) layer 3, c) layer 4, and d) faults plan view .............. 106 Figure 5.8 Transient state simulated- versus observed heads ............................................................ 107 Figure 5.9 Transient state residuals in calibration- (left) and validation phase (right) ...................... 107 Figure 5.10 Pipe water discharge changes for the scenarios with highest changes ........................... 111 Figure 5.11 Layout of hypothetical- and adjusted pipes. ................................................................... 116 Figure 5.12 Water leakage changes for hypothetical pipes with the largest effects .......................... 118 Figure 5.13 Layout of three connected pipes in three layers in the left abutment ............................ 120 XII List of Tables Table 1.1. Some Karun 4 dam technical specifications .......................................................................... 6 Table 2.1. Frequency of Lu-values within the various classes for the different borehole groups ........ 25 Table 2.2. Frequency of grout-Take within the various classes for the different borehole groups ...... 26 Table 2.3. Frequency of Lu-Take pairs within the four Ewert’s groups for the various boreholes ..... 28 Table 3.1. Karun 4 dam technical specifications .................................................................................. 35 Table 3.2. Standpipe piezometers results of the left and right banks of the grout curtain ................... 43 Table 3.3. Results of observation wells’ measurements ....................................................................... 44 Table 4.1. Lithological specifications of Asmari 1 subunits ................................................................ 54 Table 4.2. Layer specifications assigned to the aquifer model ............................................................. 61 Table 4.3. Proposed hydraulic conductivity values by Sadd Tunnel Pars company (2015) ................ 64 Table 4.4. Observation wells’ and piezometers’ specifications ........................................................... 66 Table 4.5. Stress periods with dates and lengths of transient model .................................................... 69 Table 4.6 Characteristics of the different sections of the drainage curtain. ......................................... 73 Table 4.7. Hydraulic conductivity (m/day) of the grout curtain derived from the calibrate model ..... 77 Table 4.8. Pipe specifications assigned to the CFP- model .................................................................. 82 Table 5.1. Lithological specifications of Asmari 1 subunits ................................................................ 97 Table 5.2. Layer specifications assigned to the aquifer model ........................................................... 103 Table 5.3. Finally calibrated pipe specifications (K- and D- values) of the CFP- model ................ 105 Table 5.4. Statistical results of the calibration and validation of the model ....................................... 107 Table 5.5. Cumulative volume budget in (m3) for the entire model for RS ...................................... 108 Table 5.6. Changes (%) for the different pipe diameter scenarios relative to RS ............................. 110 Table 5.7. Scenarios with the highest changes of pipes water discharge relative to RS .................... 111 Table 5.8. Scenarios with highest drainage differences, relative to that of the RS. .......................... 112 Table 5.9. Geological specifications of hypothetical pipes in the three layers (bottom to top) ......... 115 Table 5.10. Changes of water escape in comparison to RS for 3 hypothetical pipes scenarios ......... 117 Table 5.11. Results of different scenarios of the connected pipes. .................................................... 121 XIII List of Abbreviations and Acronyms AS1 Asmari 1 geological formation C Conductance of drainage curtain (m2 per day) CFP Conduit Flow Process for MODFLOW-2005 CHD Time-variant specified head package (MOFDLOW-2005) D Diameter (meter) DRN Drainage package (MODFLOW-2005) EC Electrical conductivity (here: μS/cm or micro Siemens per centimeter) HOB Head Observation package (MODFLOW-2005) IWPCo Iran Water and Power Resources development Company K Hydraulic Conductivity (meter per day) LG Left Gallery of the grout curtain LOW Left abutments Observation Well Lu Lugeon Value m.a.s.l mean water above sea level (meter) OW Observation Well P standpipe Piezometer R2 Coefficient of determination or R-Square Re Reynold’s number (dimensionless) RG Right Gallery of the grout curtain ROW Row abutment Observation Well rmsr Root mean square residual RQD Rock Quality designation RS Reference Scenario Ss Specific Storage Sy Specific Yield 1 General Introduction I.Background During dam construction and later in the operating phase, some geological, hydrogeological and geotechnical problems can take place. One of the most challenging of the latter is water leakage. Leakage not only causes water loss of the reservoir, but can, at the same time, threaten the dam stability which, eventually, may lead to dam break. One of the potentially paths for water leakage in dams is fractured and faulted zones, which could conduct water flow from upstream to downstream. Flow theory in a classical porous media was presented by Darcy more than 150 years ago based on experimental studies he did. But formulating fluid movement theory in fractured and faulted zones has been facing a lot of problems because of different size of empty spaces between particles which precludes the use of classical porous media flow theory. Classification of the karst aquifers from the porosity point of view are: - Primary porosity which is established during sedimentation and diagenesis. Karst aquifers store a small amount of water normally. - Secondary porosity is fractured and faulted zones which has been formed by tectonics forces. - Third porosity that has been formed by solution processes. This multiplicity of porosity causes complicated flow regimes in karsts. Lack of precious data of conduit geometry has intensified such complications. The Karst aquifers have been classified based on their hydraulic specifications to following categories: - Mature system with a conduit system that has a high hydraulic conductivity (for instance Alps karst system) and which have been formed by fracturing and faulting processes or by the dissolution of conduits. In such a media groundwater flow is essentially a turbulent pipe flow as the pore openings are more than 1 cm wide. - Dispersion system which has a conduit matrix with openings less than 1cm wider Groundwater tends to behave as linear flow. - Hybrid system which is a combination of two above mentioned systems. In this case both conduit and matrix systems exist. Most karst aquifers can be classified as a hybrid system, whereby the karst features lead to an aquifer with high transfer velocities (high hydraulic conductivities) and significant storage of 2 water. Therefore, any attempt of karst aquifer modeling must take into account the special porosity structure of such a media, namely, the so-called double continuum porosity. II.Case study: Karun 4 dam / Iran The Karun 4 dam is Iran highest dam. It was constructed over a karst limestone and water leakage is one of the most important problems of this dam. For more information please see: http://en.iwpco.ir/Karun4/default.aspx. The Karun-4 dam is located on the Karun river in the state of Chaharmahal and Bakhtiari, at a distance of 180 km from the State Capital of Shahrekord, Southwest of Iran. The dam is located in the tectonically folded belt of the Zagros mountains within a canyon that is located in the southwestern flank of Kuh-sefid anticline with a NW-SE trend. Most of reservoir is situated on Cretaceous to Miocene limestone and marly limestone. The geological formations around the dam site consist of carbonate layers of Asmari formation (As) which can be divided into lower (As1), middle (As2) and upper (As3) units and marlstone and marly limestone of impervious Pabdeh formation (Pd) on upstream of dam site. AS1 outcrops in dam valley abutments and are highlighted in the cliff around the site. The strike of the beds is nearly parallel to the dam axis and the dip is at an angle 70-80 degrees in downstream direction. The AS1 unit, with an apparent thickness of about 200 m in normal direction to the bedding, constitutes the foundation rock of the dam. The AS1 unit is divided to more detailed 10 units from AS1-a to AS1-j, from upstream to downstream. Bedding thickness of the rock mass varies from 2 m to 30 m. Table 1 shows subunits of AS1 and description of their lithological units. The geology of the dam site is affected by many faults, 122 major joints and 4 discontinuities sets. The joint openings along the bedding at the ground surface range from 1 to 100 mm and decrease with increasing depth. III.Research objectives The ultimate purpose of this research is the numerical simulation of the hydraulic effects of the dam foundation and abutments, to determine the paths of the water leakage at the dam site, as well as to estimate the effects of the dam impounding on the water resources downstream of the dam. The research proposed is of particular importance for the country of Iran, the more so, since the Karun-4 dam as Iran highest dam has a strategic place in the water and power 3 national network, so that any instability- or safety problems with this dam would result in irrecoverable losses to the country as a whole. IV.Structure of dissertation This PhD dissertation is written in five chapters to demonstrate the effects of construction of the Karun 4 dam on the groundwater flow in the adjacent aquifer and investigate the groundwater flow regime in a karstic dam site using a numerical model which simulate the dual porous media of the dam site. The specific research objectives of this thesis are as follows: - In chapter 1, the effects of the construction of the Karun4 dam in south-western of Iran, on the groundwater level in the adjacent aquifer is estimated, using Lugeon Test and Chemical investigations results from the dam site before and after impounding of the reservoir. - The grout curtain implemented on the dam foundation and abutments of Karun 4 dam prevents the water leakage from the reservoir. In chapter 2 the permeability of the grout curtain is evaluated using the Lugeon values and the grout take amounts. - Chapter 3 presents the seepage controlling system of Karun 4 dam which includes the grout curtain, drainage curtain and monitoring instruments like standpipe piezometers and observation wells. The monitoring results of this system indicate the need of numerical modelling. - Chapter 4 describes the simulation of groundwater flow at the dual porous media of Karun 4 dam foundation and abutments using the Conduit Flow Process (CFP) for MODFLOW-2005. The acquired data from dam site used to calibrate and validate the model. The laminar and turbulent flow regime are modelled simultaneously. - In chapter 5, the different scenarios of the water leakage through the dam foundation and abutments are investigated using the calibrated numerical model. The critical faults from geological surveys are identified, and the possible increasing water losses are discussed. 4 Chapter 1- Construction Effects of the Karun 4 Dam, Iran, on the Groundwater in the adjacent Karstic Aquifer Abstract The Karun 4 concrete dam is with a maximum dam height above foundation of 230m the highest dam of Iran. It has been in operation since July 2011. It was constructed over a karst limestone and water leakage has been, and continues to be, one of the most important problems of this dam. The geological formations around the dam site consist of carbonate layers of the Asmari formation and marlstone and marly limestone of the impervious Pabdeh formation, upstream of the dam site. The object of this study is to determine the dam construction effects on the groundwater levels in the adjacent aquifer. Thus, before the dam construction, the natural groundwater levels at both banks were recorded to be 5-8m higher than the river stage. This indicates that the natural groundwater gradient at that time was from the banks towards the river, i.e. the latter was an effluent stream before the construction of dam. At that time no important springs were identified at the dam site., other than a minor spring with a discharge <5 lit/sec that emanated at the contact interface between the permeable Asmari- and the impermeable Pabdeh formation (perching aquifer conditions) on the left bank, upstream of the dam. A Lugeon test that was carried out at that time indicated that the permeability of the adjacent limestone is high, as it varies from 25 to 55 Lu in the layers above the river level, but it decreases gradually to 3 Lu in the formations below. After dam impounding, some changes in the borehole's water levels were observed. Thus it has been found, in particular, that the leakage from the reservoir has induced groundwater level rises between 12 to 17 meters. Keywords: Karun 4 dam, Iran, water seepage, karst, hydrogeology 5 1.1 Introduction When it comes to the construction of dams and reservoirs in karsts, which are known for their high perviousness due to dissolution faults and conduits, the prevailing risk is water loss, while the risk of stability is considered to be considerably lower (Milanović, 2004). The term Karst is generally defined as a terrain underlain by rocks that are highly soluble (such as limestone) with well-developed secondary porosity, and a terrain that exhibits a distinctive hydrology and landforms such as sink holes, sinking streams, springs and caves (Ford & Williams, 2007). Water loss through basements and abutments of the dams constructed in karst areas leads to considerable costs and non-productivity of the dams. Although the construction of a grout curtain is the most common method to reduce or prevent the water loss, it can solve the problem only partly and in some cases this operation is even fruitless. For example, the Lar Dam in Iran has used less than half of its conceived capacity, because of continuous leakage - about 10 m3/s - through the dam foundations and abutments, since its construction in 1980 (Uromeihy, 2000). Although efforts, such as water injection into the fault, lining of parts of the bottom of the reservoir by thick blankets of clay, or the construction of cut off walls, have been made in recent times to prevent or, at least, to reduce the water leakage, successes to that regard have been very limited. The main paths of water loss underneath dams and reservoirs in karst are the dissolution conduits and the main cause of this problem then is the poor identification of such conduits as being responsible for the leakage. The generation and development of karst conduits within a rock mass is not completely random, but is somehow guided by various chemical pre-requisites which drive the dissolution processes (Wright, 1991). In fact, the primary flow paths in karst rocks are discontinuities such as joints, faults and bedding planes (Kiraly, 1975), although intrinsic porosity plays also a non-negligible role ((Bakalowicz, 2005). The Karun 4 concrete dam, which is the focus of the present study, is with a maximum dam height above foundation of 230m the highest dam of Iran. This dam sits on a carbonate (limestone) rock formation, so the problems mentioned above to exist in such a geological environment are to be expected. After the completion of the dam, the Karun 4 reservoir reached the full water impound phase in 2010. The object of this study is to investigate the likely effects of the reservoir’s water impound on the groundwater underneath the dam site. More specifically, this paper is concentrating on results of the analyses of various hydrogeological and hydrochemical investigations and, in particular, of Lugeon tests. 6 1.2 Study area 1.2.1 Site specification and geology Karun is the largest river and the only navigable waterway in Iran. It originates in the Zagros mountains, west of Isfahan, flows out of the central Zagros range, traverses the Khuzestan plain, and joins the Shatt al-Arab, with the latter then discharging into the Persian Gulf. There are a number of dams on the Karun river, such as the Gotvand-, Masjed Soleyman-, Karun-1 (Shahid Abbaspour)-, Karun-3 and Karun-4 dam, all mainly built to generate hydroelectric power and to provide flood control. The Karun 4 Dam itself is located on the Karun River in the state of Chaharmahal and Bakhtiari, at a distance of 180 km from the State Capital of Shahrekord, in southwest Iran (Figure 1.1). The objectives of the Karun 4 Dam construction have been the following: - regulation of the Karun river discharge of 3.7 billion m³/year - control of destructive floods and overflows of Karun River - generation of hydroelectric energy of about 2107 billion kWh per year. The dam has been in operation since July 2011, after termination of the water impound phase, which started in 2010. Table 1.1 shows some technical specifications of the dam. Table 1.1. Some Karun 4 dam technical specifications Dam type Concrete double arc Height from foundation 230 m Crest length /width 440 m / 7m Foundation width 37 to 52 m Total volume of concrete placement 1.65 million m3 Total volume of the reservoir 2190 million m3 Power Plant Type Ground Spillway Chute, orifice types and free The dam is located in the tectonically folded belt of the Zagros Mountains within a canyon that is located in the southwestern flank of the Kuh-Sefid anticline which has a NW-SE trend. Most of the reservoir is situated on Cretaceous to Miocene limestone and marly limestone (Haghnejad et al., 2014). The geological formations around the dam site consist of carbonate layers of the Asmari formation (As) which can be divided into lower (As1), middle (As2) and upper (As3) units and marlstone and marly limestone of impervious Pabdeh formation (Pd) on 7 Figure 1.1. Location of Karun 4 dam in Iran upstream of dam site (see Figure 1.2 and Figure 1.3). The lower unit (As1) comprises limestone, porous limestone and some marly limestone with some marlstone, generally thick to very thickly bedded and highly karstic. The measured mean RQD (rock-quality designation) index (Deere, 1964) ranges between 55% and 83%. The middle unit (As2) is located downstream of the dam and comprises alternating limestone, dolomitic, marly limestone and marl. Karstification is limited and the RQD is between 53% and 84%. The upper part (As3) consists of alternating layers of limestone, marly limestone and marl with slight karstification and an RQD between 45% and 78%. These RQD-values indicate a fair to good quality rock in all three units (Koleini, 2012). As1 outcrops in the dam valley abutments are highlighted in the cliff around the site. The strike of the beddings is northwest to southeast, and southern limb slopes to the south with a dip of 70-80 degrees. The AS1 unit, with an apparent thickness of about 200 m in normal direction to the bedding, constitutes the foundation rock of the dam. The geology of the dam site is affected by many faults, 122 major joints and 4 discontinuities sets. The strongly dipping reverse Monj fault, with a length of 14 km, is located at a distance of about 500 m downstream of the dam. The joint openings along the bedding at the ground surface range from 1 to 100 mm and decrease with increasing depth. 8 Figure 1.2. Karun 4 dam geological plan Figure 1.3. Karun 4 dam geological cross-section viewed from the reservoir downstream towards the dam 9 1.2.2 Hydrogeology Before dam construction in year 2002, the natural groundwater level at both banks was recorded to be 5-8m higher than that of the river water level. This indicates that the natural gradient of groundwater at that time was from the banks toward river, i.e. the river was an effluent stream, which is the normal situation for most rivers in the world. In fact, the phreatic levels in the left and right banks of the river were measured as 848 and 845m NN, respectively. In fact, at that time, before 2002, no important springs were identified at the dam site. The only minor spring, with a discharge of less than 5 lit/sec, emanated at the contact interface between the permeable Asmari Formation and the impermeable Pabdeh formation, upstream of the left abutment of the dam. In the Asmari Formation dissolutional surface features such as vugs and small cavities along the discontinuities with infilling were observed. Some larger cavities (1 to 2m in diameter) were encountered during the excavation of the galleries. A porous limestone layer was identified, which had a vuggy feature and was recognized as a key bed underneath the dam site. Lugeon tests were carried out at that time (2002), which resulted in a rather high permeability between 25 to 55 Lu in the limestone layers above the river level, but which decrease gradually to 3 Lu, as one goes towards the formations below the river water. After dam impounding in 2010, some changes in the borehole water levels were observed and it was found that the leakage probability from the grout curtain had increased. This was corroborated by the studies of Shabab (2011) and Kamali and Ashjari (2012) who proved the existence of increased leakage, but these authors could neither specify nor identify the pathways of the leaking waters. 1.3 Study method Geology maps of the dam site were prepared, based on geology maps of 1/100000 Geology Surveys of Iran and 1/20000 of Mahab Ghodes Consulting Company. The geotechnical data of drilled boreholes, such as geological logs, RQD, Lugeon tests, and placement were collected. From the Water Resource Management Office of Iran data about the following hydraulic and hydrogeological features were gathered: the flow system and hydrogeological characteristics of the sites, water level measurements in boreholes in karstic rock mass or wells, locations of springs and their discharges. With the geological logs and geotechnical data potential horizons of conduit developments were found. The Lugeon tests results identified high permeability zones in a borehole that could hint of leakage zones in the geological formations under the 10 dam. Finally, temporal variations of ground water levels and several chemical parameters were correlated with the height/volume of the water impound in the reservoir. 1.4 Results 1.4.1 Before water impound Lugeon tests and a hydrogeological assessment of the area were done before the dam water impounding in 2010, in order to evaluate the permeability properties of the karst formations at the dam site (Mahab Ghods report, 2002) and get some hints of possible leakage pathways. The results of these tests which were performed in boreholes drilled at a ground surface height, ranging between 867 and 885 m NN, indicate a high degree of heterogeneity in both vertical and horizontal directions in the left and right abutments and the foundation of dam. In the left abutment, the permeability variations define three major zones: 1) The carbonate rocks containing marl. These rocks have fractures and are very porous, due to faults and joints which cause a strong vertical variation of the permeability. A high value of the latter of more than 100 Lu, has been reported in a very deep seated layer, at depths between 193 and 200 m. But the thickest highly permeable zone lies in the depth range of 60 to 120 m. 2) Places with high variations of the permeability in a shallow depth range, but a consistently thick layer of very high permeability between 135 to 160 m depth. 3) Boreholes bottoming near to the contact zones of the named Asmari and Pabdeh formations show the lowest permeability in the left abutment of the dam body, with values lower than 3 Lu in the Pabdeh part of boreholes. In the right abutment of the dam body, the permeability variation is more consistent with the lithology of the layers than in the left one. Thus, the permeability is very low, less than 8 Lu, in a depth of 100 m and below, whereas high permeability zones are located in depths less than 60 m. The permeability of rocks beneath the foundation itself is very low, with the highest Lugeon values of only 8 to 32 determined in the depth range of 10.5 to 44 m , and of even less than 8 in other horizons. The water levels of the Karun and Monj rivers are 806 and 805 above NN, respectively. The groundwater levels of karst aquifer varied from 852 to 855.5 at the dam site in 2002 (Figure 1.3.), which means that they were higher than the Karun river stage and lower than that of the 11 Monj river. Hydrogeologically, this confirms that at that time the Karun river acted as a gaining- and the Monj River as a loosing stream. 1.4.2 During and after water impound The reservoir water impounding started on March 25, 2010. Up to date it has not yet reached the target level, due to a series of problems related to the dam structure and the named water seepage. The following analysis is limited to data of year 2010, because of the inaccessibility of more recent data. During the six months after the beginning of the impounding phase (April – October 2010) the reservoir water levels as well as groundwater levels in several boreholes of the grout and the drainage galleries were recorded on a weekly interval. Figure 1.4. shows these reservoir water levels together with groundwater levels in two selected boreholes of the left and right abutments. From April- October 2010, the reservoir level increased steadfastly from 806 to 920m NN. In fact, a final recording at the end of that year provided a value of 946 m NN for the water level in the reservoir (Kamali and Ashjari, 2012). One may notice that during the April- October 2010 time period the groundwater levels in the boreholes of the left and right abutments show some different behavior. Thus, whereas in the right abutment a hydraulic head increases of up to 20 m is observed in the first three months, leveling off thereafter, the situation is opposite in the left abutment, as the groundwater level starts to rise only in the second half of the 6 months observation period. In general, the groundwater level variations in the right dam abutment turned out to be more complicated and depend on the distance of the borehole from the junction of the dam with the drainage galleries. Thus, at a distance of less than 50 m, the head changes are less than 9 m, going up to 20-24 m in the distance range 50-144 m, and, eventually, reaching 15-17 m in the 44-250 m range. Therefore, the head variations in the right abutment are more responsive to the reservoir water impounding than those in the left one. As for the boreholes in the abutment close to the dam (not shown here), only low groundwater head increases of about 5 m were measured during the named time period. On the other hand, at distances of 114 to 138 m upstream of the junction of the dam with one of the galleries in the left abutment, the increases are with 20 to 30 m the highest, but they become lower again further away. 12 Figure 1.4. Change of water levels after water impounding of Karun 4 dam (2010) 1.4.3 Hydrochemistry of the groundwater Since 2010, after commencement of the water impounding, electrical conductivity (EC) and major ions including CO3, Cl, Na, Ca, K, Mg and SO4 have been measured on a monthly basis in the boreholes as well as in downstream sections of the Karun and Monj rivers. The temporal variations of these parameters over a time period of nearly four years are shown in the various panels of Figure 1.5. Regarding the EC, one can notice from the corresponding plot that the EC of the Karun river has been varying in a rather narrow range of 400 to 560 µS/cm since the time of water impounding and it has been consistently lower than that of the Monji river, up to the end of 2011. After that time the EC- values and trends in both rivers are similar, owing to the fact that the impounded reservoir started to extend to the Monj river catchment. The ECs of the left abutment are high in the beginning of the water impounding phase and decrease slowly with the progressive water level increase in the reservoir. They become then closer to the EC-values of the Karun river. It should be noted that the first few records of the EC in the right abutment are missing, so that the EC-variations at the beginning of the water impounding process are not known. Nevertheless, it appears that the EC- decreases in the right abutment are similar to those of the left one . However, because of the higher permeability in the right than in the left abutment, it 13 is of no surprise that the EC in the former decreases also faster than that of the latter towards the value of the river water. The hydrochemistry of the groundwater and two rivers around the Karun dam site consists of mainly bicarbonate or bicarbonate-sulfate ions, owing to the contact of these waters with the named Asmari and Gachsaran limestone formations. Now, if one assumes the possibility of leakage of reservoir water into the aquifer, the chemical properties of the water in the reservoir boreholes and the reservoir should become similar after some time. In fact, the concentration plots of the various ions in Figure 1.5. show, overall, a decreasing trend of the borehole concentrations in the beginning of the water impounding phase, after which they follow fairly well the trends of the corresponding parameters measured in the Karun River. Therefore, these chemical variations of boreholes provide evidence for water seepage from the dam site. These trends are weak for the bicarbonate (CO3) variations, due to the effect of carbonate formation and dissolution along the water pathways though the limestone layers. 1.5 Conclusions and outlook The following results obtained from the various hydrological and hydro-chemical analyses show that water seepage from the dam reservoir to the aquifer has occurred since the commencement of water impounding in the Karun 4 dam reservoir in early spring of 2010: 1) Groundwater levels have increased significantly by up to 20 m in both abutments of the dam. 2) The left abutment shows a lower hydraulic head increase than the right one. 3) With ongoing water impounding, the concentrations of the major ions and as well the EC of the groundwater get closer and closer to those in the reservoir. The major reason of such a response of the groundwater to the water impounding in the reservoir is the high karstification of the limestone at the dam site, as indicated by the following factors: 1) Lithologically, the dam site consists of pure limestone which is highly prone to dissolution. 2) Caves, which are a main surface indicator of karstification in a region, have also been reported at the Karun 4 dam site. 3) The faults and fractures created during the karstification process provide suitable water pathways from the dam site to the aquifer. 14 Figure 1.5. Variations of hydrochemical parameters between 2010 and 2014 15 4) Lugeon test results show high permeability values in different horizons of the karst terrain around and underneath the dam site. Today (2014) the dam is going to be repaired, due to improper grout curtain installation in the abutment and the foundation. In fact, based on the various indices of karstification in the region, the mentioned problems with the Karun 4 dam were already partly predicated in some of the early planning reports (Kamali and Ashjari, 2012), however, the recommendations of various geologists were not adhered to, with, in addition to an unsuitable monitoring of the grout curtain installation, has led to the series of problems encountered at the dam at present time. The next step in this ongoing study is the development of a groundwater flow and transport model, capable of simulating water leakage and solute transport through the karstic formations at the dam site. Much research along these lines is still needed, as the numerical flow and transport models for non-Darcian porous media like in karst, which consists of a mixture of pores and conduits are still lacking reliability (Bakalowicz, 2005). Various model approaches have been proposed to tackle this problem, the most notable ones known under the notation of double continuum (porosity) models, which represent the fractured limestone as a mixture of a conduit matrix and a classical porous medium, with appropriate exchange terms (Teutsch and Sauter, 1998; Scanlon et al., 2003). The outcome of the research proposed further is of particular importance for the country of Iran, the more so, since the Karun 4 dam, as Iran’s highest dam, has a strategic place in its water and power national network, so that any instability- or safety problems with this dam would result in irrecoverable losses to the country as a whole. References Bakalowicz M (2005). Karst groundwater: a challenge for new resources, Hydrogeol J., 13,148–160 Deere DU (1964). Technical description of rock cores, Rock Mechanics Engineering Geology, 1, 16- 22. Ford D & Williams P (2007). Karst Hydrogeology and Geomorphology. Wiley Ltd, England. Haghnejad A, Ahangari K & Noorzad A (2014). Investigation on various relations between uniaxial compressive strength, elasticity and deformation modulus of Asmari formation in Iran. Arab. Journ. Science and Engineer., 39, 4, 2677-2682. 16 Kamali M & Ashjari J (2012). Determination of efficiency of grout curtain of Karun 4’s Dam by dye tracing test. Report to IWPC (Iran Water & Power resources development Company), 124 p. Kiraly L (1975). Rapport sur l'état actual des connaissances dans le domaine des caractères physiques des roches karstiques, in: Bueger and Dubertret (Eds.), International Union Geol. Sci., Series B. 3, pp.53-67. Koleini M (2012). Engineering geological assessment and rock mass characterization of the Asmari formation (Zagros range) as large dam foundation rocks in southwestern Iran. Ph.D thesis, University of Pretoria, South Africa. Mahab Ghods Consulting Engineering Co. (2002). Karun 4 dam project. Final engineering geological report (phase 2). Milanović PT (2004). Water resources engineering in karst. CRC Press, Boca Raton. Scanlon BR, Mace RE, Barrett ME & Smith B (2003). Can we simulate regional groundwater flow in a karst system using equivalent porous media models? Case study, Barton Springs Edwards aquifer, USA. J. Hydrol., 276, 137–158. Shabab Boroujeni B (2011). Determination of effective factors on permeability in karstic dams in west of Zagros. M.Sc. thesis, Tehran University, Iran. Teutsch G & Sauter M (1998). Distributed parameter modelling approaches in karst-hydrological investigations. Bull. d’Hydrogeologie, 16, 99–109. Urumeihy A (2000). The Lar Dam; an example of infrastructural development in a geologically active karstic region. Journal of Asian Earth Sciences. 18, 25-31. Wright VP (1991). Palaeokarst types, recognition,controls and associations, in: Wright, V.P., Esteban, M., Smart, P.L. (Eds.), Palaeokarsts and Palaeokarstic Reservoirs, Reading: P.R.I.S. Occasional Publication Series 2, pp.56-88. 17 Chapter 2- Evaluating Permeability and Groutability at the Karun 4 Dam, Iran, using Lugeon Values and Grout Take Abstract Karun 4 dam is the highest concrete double arc dam in Iran located in the Zagros mountain range in southwest of Iran. It has been under operation since March 2011. The dam foundations and abutments are within the Asmari formation, consisting of fractured carbonate rocks with a high potential of karstification. The water tightness of the dam was one of the main issues during the construction. The grout curtain was constructed to solve the problem in three steps. The object of this study is to evaluate the designed curtain based on the grout taking rate, permeability test, and geological logs. The studied boreholes consist firstly of the exploration boreholes, which represent the structural situation and characteristics of the rock formation, and, secondly of the check holes that were drilled after finalizing the implementation of the grout curtain. The results illustrate the efficiency of the grout curtain, but, nevertheless, some water leakage is still occurring through the curtain, so that further monitoring and treatment is required to ensure the dam safety, namely, to prevent the water leakage. Keywords: Permeability, Groutability, Lugeon test, Grout take, Water leakage, Grout curtain, Karst, Karun 4 dam, Carbonate aquifer. 18 2.1 Introduction Water leakage through the foundation and abutments of a dam which is built on limestone formations, is one of the most important challenges of the big dams’ construction and operation. The water- soluble carbonate structure of such limestone may lead to the development of fissures and fractures that may, eventually, expand to conduits and even caverns and caves; a process which is known as karstification (Milanović, 2004; Ford & Williams, 2007). As a consequence, water from a dam reservoir can easily seep through its foundation and abutments and so ultimately affect the dam’s stability. Extreme examples for such a situation to occur are the Camaras dam in Spain (Jaeger 1979) and the Lar dam in northern Iran (Uromiyehi 2000). One of the most prevalent methods to reduce or even to avoid this adverse process of water leakage is to construct a virtual curtain by drilling a very intensive borehole-spreadsheet in the dam foundation and abutments, and then injecting cement grout (Fell et al. 2015). To evaluate the function and efficiency of such a grout curtain, the permeability or hydraulic conductivity of the mass rock should be measured before and after its construction. A test to do this has been proposed first by Lugeon (1933), known as the Lugeon test, is still up-to-date the most widely used one for that purpose. However, for evaluating the permeability of a rock mass in more objective manner, it has been strongly recommended to take into account, in addition to the Lugeon values, also the geological and lithological characteristics, cement Take, etc. (Goodman et al. 1964, Barton et al. 1985). Ewert (1985) analyzed in detail the Lugeon test used in rock grouting, with emphasis on dam sites. Houlsby (1990) published the “Construction and Design of Cement Grouting” reference book and classified the relationship between the amount of cement -Take (from now onwards: Take) and the hydraulic conductivity inferred from a Lugeon test into four categories. Foyo et al. (1990) modified the water pressure test to analyze the permeability characteristics and later Foye et al. (1997) studied geological features, permeability and groutability-characteristics of the Zimapan Dam foundation in Mexico. Using a Back-Propagation Neutral Network (BPN), Yang (2004) analyzed the relation between Lugeon test values and cement Take in grout holes of curtains of the Li-Yu-Tan dam in Taiwan. Foyo et al. (2005) proposed a Secondary Permeability Index (SPI), obtained from water pressure tests in dam foundations. They claimed that SPI is more quantitative than classical Lugeon values. Uromeihy und Barzegari (2007) studied the seepage problem of Chapar-Abad dam in Iran with Lugeon values and seepage 19 modelling using the PLAXIS software. Bonacci and Roje-Bonacci (2008) studied the water losses of the Ricice reservoir built in the Dinaric karst in Croatia and Sadeghiyeh et al. (2013) compared the permeability (Lugeon) and groutability (cement Take) at the Ostur dam site rock mass in Iran to design the appropriate grout curtain. Azimian and Ajalloeian (2013), referencing Foyo et al. (2004) studies, indicated that SPI and Lugeon describe more or less the same phenomenon, and differ only by a proportional constant which depends on the borehole diameters. Similar studies were carried out by Kocbay and Kilic (2006) at the Obruck dam in Turkey; Ghafoori et al. (2011) at the Daroongar dam in Iran, Gurocak and Alemdag (2012) at the Atasu dam in Turkey; Uromeihy and Farrokhi (2012) at the Kamal-Saleh dam in Iran; Berhane and Walraevens (2013) at the Geba dam in Ethiopia and Rajabi et al. (2015) at the Tanguyeh dam in Iran. In the present paper, the permeability as the most important parameter determining water seepage of the Karun 4 dam - constructed on a limestone formation with karst potential (Mahab Ghods report 2002) – will be studied. More particularly, experimental results from various engineering geological investigations at the dam site, such as measured water pressure tests and cement Take are used to determine the permeability in numerous boreholes and check holes across the dam foundation and abutments. The focus is on a comparison of the groutability and permeability in various holes of the grouting galleries using Ewert’s (1985) proposed method. As mentioned above, the particular highlight of the present study is that, by being able to compare the measured permeabilities before and after construction of the grout curtain, a more comprehensive view of the efficiency of the grout curtain and, consequently, of the proper function of the dam reservoir can be obtained. 2.2 Study area and dam specifications 2.2.1 Site specification and geology The Karun 4 dam is located close to the source of the Karun River in the state of Chaharmahal and Bakhtiari, in southwest Iran (Figure 2.1) and has been impounded since March 2010. Its purpose is to regulate the Karun river’s discharge and to control the frequent destructive floods of the Karun River, but also to generate about 1GW hydroelectric energy. 20 Figure 2.1. Location and details of the Karun 4 dam in Iran The Karun 4 dam is in the tectonically folded belt of the Zagros Mountains within a canyon that is located in the southwestern flank of the Kuh-Sefid anticline which has a NW-SE trend. Most of the reservoir is situated on Cretaceous to Miocene limestone and marly limestone (Haghnejad et al., 2014). The geological formations around the dam site consist of carbonate layers of the Asmari formation (As) which can be divided into lower (As1), middle (As2) and upper (As3) units and marlstone and marly limestone of impervious Pabdeh formation (Pd) upstream of the dam site (see Figure 2.1). The geology of the dam site is affected by 11 faults in the left abutments and 9 faults in the right abutment. The joint openings along the bedding at the ground surface range from 1 to 100 mm and decrease with increasing depth. The lower unit of the Asmari formation AS1 consisting of limestone, marly limestone and marlstone with some porous limestone (Hosseiny Sohi et al, 2014), subdivided from upstream to downstream into 10 more detailed subunits AS1-a to AS1-j. Furthermore, a porous limestone layer (AS1- g) was identified, which had a vuggy feature and was recognized as a key bed underneath the dam site. 2.2.2 Dam specifications and grout curtain To seal the dam foundation and abutments, a very large grout curtain was designed and constructed. The construction was done through the excavation of 5 series of galleries in both 21 Figure 2.2. Selected Exploratory- , line 1 and line 2 check holes in the grouting galleries, View from reservoir toward downstream. To name the boreholes, the following notations are used, for example, borehole L1A024, left gallery 1, LG1, exploratory hole (A), at a distance of 24 m from the entrance of the gallery, and R4CH248 for a checkhole in the right gallery 4 at a distance fo 248m, and 2R4CH114, for a check hole of the grout curtain (marked with number 2) left and right abutments at different elevations from the dam crest elevation to the bottom and, finally, the foundation gallery, called the 806-gallery (see Figure 2.2). From inside these galleries the grouting boreholes are drilled, examined and injected, to have an integrated virtually wall or barrier, called grout curtain against water leakage. Some of the primary drilled holes were selected as exploratory boreholes and the rock permeability tests and geological logs were carried out within these. They were then considered in the grouting phase as normal drilled boreholes and the amounts of their cement Take were reported. Based on the obtained results from these exploratory boreholes, the grout curtain in two lines was designed and constructed, to sew/stitch the dam structure to the Pabdeh impermeable formation in the upstream. The first line was done almost along the galleries, whereas the second line, located 0.9 meter upstream toward the reservoir, was drilled principally in the early sections of the galleries, because after grouting the first line, the farther sections of the galleries exhibited very few permeability and there was thus no need to grout them in addition. After implementation of the first and second lines of the grout curtain, the line 1 and 2 check holes were drilled, and the permeability and Take were measured. The seepage controlling system consists of the grout curtain and the drainage curtain. The latter consists of boreholes drilled at 40 meters’ horizontal distance toward downstream in two galleries in the left and right abutments and allows not only to detect and to explore the still leaking water through the grout curtain, but also aids to reduce the uplift pressure and so prevents the breaking of mass rock downstream. 22 2.3 Methodology and data 2.3.1 Water pressure test (WPT) or Lugeon test Water pressure test, packer test or simply Lugeon test has been developed by Professor Maurice Lugeon (1933). A Lugeon unit is defined as one liter/minute of water absorption per meter of test length of drill hole when the water in the borehole remains at a pressure of 10 bars or 1 MPa over a period of 10 minutes (Lugeon, 1933; Houlsby, 1990; Singh and Goel, 1999). 2.3.2 Take The (grout) Take is defined as the rate of dry cement mass that is injected to the test section, divided to the length of section (kg/m) (Deere 1982). 2.3.3 Ewert method Ewert (1985) proposed a qualitative method to interpret the relationship between permeability and groutability in the dam foundations. He classified the obtained results from different dam sites and tests in to Figure 2.3 grouping: Group A: large amount of Lugeon value and low grout Take indicates that water can pass through the numerous fine fissures of the rock, but the corresponding grout is not permitted. Group B: Approximately proportionality between water and Take justifies well the necessity of grouting. Group C: The small water pass versus large Take cause to the hydraulic fracture in the rock mass and shows that either the injection pressure is not proper or at all unnecessary. Group D: Low water and Take signs the sealing of area and no need no treatment. Comparing the four groups with each other shows that grouting is only in the Group B section necessary and reasonable. Due to small Take, a real grout curtain cannot be expected in group A, whereas in Group C and especially group D the low water leakage aim is already achieved. However, as will be discussed in the results section, due to anisotropy of the hydraulic paths and the geological properties of most rocks, namely, in limestone formations, as is the case for the Karun 4 dam, a clearly defined linear proportionality between the Lugeon values and the Take, as postulated in the Ewert method, is usually not obtained. This problem has been discussed further by Lugeon (1933), Nonveiller (1989), Kutzner (1996) and Foyo et al. (1997). 2.3.4 Data sources Mahab Ghods Consulting Engineering Company provided geology maps, layout, view and details of the grout curtain maps, the engineering geology report “Final excavation and 23 Figure 2.3. Ewert graph of water pressure test values (Lugeon) over Take injection report of the Karun 4 dam grout curtain”, rock mechanics report and the geotechnical and geological data of the drilled holes including the water pressure test data (Lugeon values) and the amounts of Take for the selected exploratory boreholes and the first and second lines of the grout curtain. In the left abutment 12 boreholes within the 5 grouting galleries at the different elevations and with different drilling depths and inclinations, are selected as representative exploratory boreholes. In the right abutment, similarly, data 16 boreholes were analysed (see Figure 2.2). Lugeon tests carried out in each borehole in any 5-meter length and then, the holes were grouted and respective Takes were measured. For line 1, data from 14 check holes in the left bank and 10 check holes in the right bank; and for line 2, 7 check holes in the left, and 5 check holes in the right bank of the grouting curtain are used (see Figure 2.2). 2.4 Results 2.4.1 Lugeon permeability The Luegon (Lu) values of the water pressure test results of the exploratory line 1- and line 2- check holes of the grout curtain are grouped separately for the left and right abutment of the dam in Table 2.1, following the classification of Ghafoori et al. (2011). Moreover, the Figure 2.4 shows corresponding barplots of the frequency distribution of the Luegon values. 24 One can notice immediately from figure 2.4 the efficiency of the grouting curtain, in as much the frequencies in the lower two Lugeon classes (0 < Lu < 3) and (3 < Lu < 10) are significantly increased in the exploratory holes, line 1- and line 2- check holes of the grout curtain for both abutments. On the other hand, only few Lu-data are within the high and very high permeability classes for the left and right abutments. These results confirm also that the permeability in the right bank is general lower than in the left bank. Based on the line 1- results, the line 2 check holes of the grout curtain were drilled in areas within the galleries where large values of Lu, i.e. strong water leakage potential, had been recorded. These positions are commonly at the beginning of the galleries from the contact point of the dam body’s concrete toward the abutment rocks, since the layout of the grouting galleries are designed to be pasted to the Pabdeh formation at the upstream of the dam. Thus, it could be expected that after implementation of the line 2 check hole of the grouting curtain, lower Lu- permeability values would be obtained. Although this is confirmed generally, in the right bank, the percentages of values of Lu in the two low-permeability classes (Lu<10) of line 1 decreases, in comparison with those of line 2. This could be explained by the locations of the drilled check holes of line 1 which are distributed along the galleries which are within the porous limestone of the different subunits of the Asmari (AS1) formation and partly within the more impermeable Pabdeh formation of, with lower Lugeon vales, whereas the check holes of line 2 extend only within the Asmari (AS1) formation. After carrying out the line 2 of grouting on the left bank, there are still 5 measurements of different check holes that reach Lugeon values more than 10 (medium permeability). They are in the porous limestone subunit of AS1-g, as the most permeable layer of the Asmari formation. The faults F4, F7 and F8 are, in the same way, potential conduits for water seepage. In fact, a total of 8 out of 10 Lu-test holes in the left bank are located either in AS1-g or close to the existing faults, which explains well the remaining medium permeability class in these places. On the right bank, two medium permeability cases are within AS1-g, with two test holes even close to the faults F15 and F13, so that the same pattern of permeability is observed here. 25 Figure 2.4. Frequency distribution of Take in the different classes for exploratory, line 1 and line 2 check holes for the left and right banks. Table 2.1. Frequency of Lu-values within the various classes for the different borehole groups 2.4.2 Groutability Following the classification of Deere (1982), Table 2.2 lists the Take results for the exploratory-, line 1- and line 2- check holes of the grout curtain and the Figure 5 shows the corresponding barplots for the left and right banks. The Take values in the exploratory holes do represent, in fact, the general groutability of the original rock mass. The positive impact of grouting line 1 can be clearly noticed from the corresponding column in table 2.2 and Figure 2.5 for all the first three groups of up to the moderately low class (Take < 50 kg/m) in both the left and right bank, whereas the values for the group of moderately high (Take>100 kg/m) and higher, decline in line 1 and line 2 check holes in both banks. Likewise, to the Lugeon tests, the Take measurements were carried out in the same boreholes. The criteria to terminate the grouting operation is to have 90% of the Take values within the - Num. % Num. % Num. % Num. % Num. % Num. % 0 - 3 Impermeable 31 20.3 52 28.1 38 33.6 53 56.4 15 27.3 19 47.5 3 - 10 Low permeability 32 20.9 62 33.5 49 43.4 32 34.0 30 54.5 15 37.5 10 - 30 Medium permeability 51 33.3 37 20.0 23 20.4 8 8.5 10 18.2 4 10.0 30 - 60 High permeability 28 18.3 24 13.0 2 1.8 1 1.1 0 0.0 2 5.0 > 60 Very high permeability 11 7.2 10 5.4 1 0.9 0 0.0 0 0.0 0 0.0 153 100 185 100 113 100 94 100 55 100 40 100Total Line 2 check holes Left Right Left Right Left Right Lugeon value Classification (Ghafoori et al., 2011) Exploratory holes Line 1 check holes 20 21 33 18 7 34 43 20 2 1 27 55 18 0 0 0 10 20 30 40 50 60 0 - 3 3 - 10 10 - 30 30 - 60 > 60 F re q u e n c y % Lugeon Values of Left Bank Exploratory holes Line 1 check holes Line 2 check holes 34 41 24 16 7 56 34 9 1 0 48 38 10 5 0 0 10 20 30 40 50 60 0 - 3 3 - 10 10 - 30 30 - 60 > 60 F re q u e n c y % Lugeon Values of Right Bank Exploratory holes Line 1 check holes Line 2 check holes 26 Figure 2.5. Frequency distribution of Take in the different classes for exploratory, line 1 and line 2 check holes for the left and right banks. <50 kg/m- class (Mahab Ghods Report, 2010). Table 2.2 shows that this goal had not yet been achieved for the check holes of line 1. In fact, that is the reason why the line 2 grout curtain was installed. And, indeed, the values in the line 2 check holes indicate that the grouting operation has been satisfactory, as now more than 90% of the Take- values are lying within the first three low grout-Take groups (see Table 2.2). Table 2.2. Frequency of grout-Take within the various classes for the different borehole groups 2.4.3 Correlation of Lu- permeability and Take- groutability based on Ewert’s method Using the method proposed by Ewert (1985) and described in Section 2.3, the Lugeon- permeability values are plotted in Figure 2.6 over the Take- groutability values for all exploratory- (top row panels) and check holes of the two grouting lines (middle- and bottom row panels) in the left (left column panels) and right (right column panels) abutments. (kg/m) Num. % Num. % Num. % Num. % Num. % Num. % 0 - 12.5 Very low 56 35,2 97 48,0 80 60,2 57 56,4 45 76,3 35 81,4 12.5 - 25 low 31 19,5 29 14,4 28 21,1 13 12,9 3 5,1 3 7,0 25 - 50 Moderately low 13 8,2 11 5,4 14 10,5 11 10,9 9 15,3 1 2,3 50 - 100 Moderate 12 7,5 14 6,9 4 3,0 12 11,9 0 0,0 1 2,3 100 - 200 Moderately high 8 5,0 14 6,9 3 2,3 4 4,0 2 3,4 1 2,3 200 - 400 High 11 6,9 21 10,4 3 2,3 2 2,0 0 0,0 1 2,3 > 400 Very high 28 17,6 16 7,9 1 0,8 2 2,0 0 0,0 1 2,3 Total 159 100 202 100 133 100 101 100 59 100 43 100 Right Left Right Left Right Left Grout Take Classification (Deere, 1982) Exploratory holes Line 1 check holes Line 2 check holes 42 17 7 7 6 9 12 59 18 11 7 3 2 1 78 6 10 1 3 1 1 0 10 20 30 40 50 60 70 80 90 0 - 12.5 12.5 - 25 25 - 50 50 - 100 100 - 200 200 - 400 > 400 F re q u e n c y % Take Values of Left Bank Exploratory holes Line 1 check holes Line 2 check holes 48 14 5 7 7 10 8 56 13 11 12 4 2 2 81 7 2 2 2 2 2 0 10 20 30 40 50 60 70 80 90 0 - 12.5 12.5 - 25 25 - 50 50 - 100 100 - 200 200 - 400 > 400 F re q u e n c y % Take Values of Right Bank Exploratory holes Line 1 check holes Line 2 check holes 27 Figure 2.6. Lugeon versus Take with the four Ewert’s groupings (see Fig. 3 for notations) for the exploratory holes (top row), line 1 check holes (middle row) and line 2 check holes (bottom row) for left (a,c,e) and right (b, d, f) abutment. 28 Figure 2.7. Frequency distribution of Lu-Take pairs in the different groups for exploratory, and line 1 and line 2 check holes for the left and right banks. Table 2.3. Frequency of Lu-Take pairs within the four Ewert’s groups for the various boreholes Following the qualitative separation of the different Ewert’s classes as indicated in Figure 2.3, the lines dividing the four Ewert- groups are also graphed in Figure 2.6. Furthermore, Table 2.3 and Figure 2.7 lists/showsthe frequency distribution of Lu-Take pairs for the various check holes in the four Ewert’s groups. Following the classifications of Lu- values and Take- values, as noted in Table 2.1 and Table 2.2, respectively, Lu- values > 60 Lu and Take- values >400 kg/m are considered to belong to the class of very high permeability. Therefore, the axis limits in Figure 2.6 have been set to these values. Figures 2.6 and 2.7 show that for the exploratory boreholes, those of the left bank have most of the Lu-Take values lying within group A which, based on the explanations of Figure 2.3, means that here high water Lu-permeability values are due to fine fissures that cannot be well grouted by more viscous grout. For the right bank, on the other hand, most of the Lu-Take values belong to group D, so that here low permeability, due to an absence of fissures, makes additional grouting unnecessary, although the vast scattering of the data here appears to warrant a further grouting treatment. Also, from Figure 2.7 and Table 2.3, it can be deduced again that, since the Lu-Take frequency in group A in the left bank is generally higher than that in the right bank, Lu- permeabilities in the former are also higher than in the latter. Num. % Num. % Num. % Num. % Num. % Num. % A 55 36,4 46 24,9 14 12,4 5 5,3 8 14,5 6 15,0 B 18 11,9 22 11,9 5 4,4 12 12,8 1 1,8 0 0,0 C 28 18,5 18 9,7 3 2,7 2 2,1 0 0,0 1 2,5 D 50 33,1 99 53,5 91 80,5 75 79,8 46 83,6 33 82,5 Total 151 100,0 185 100,0 113 100,0 94 100,0 55 100,0 40 100,0 Left bank Right bank Left bank Right bank Left bank Right bank Exploratory holes Line 1 check holes Line 2 check holes Ewert Groups 29 After grouting, in contrast, most of the Lu-Take pairs in the check holes of line 1 are shifted to the extreme (positive) group D, and this holds for both banks which shows again the good to excellent sealing capacity of this grout curtain, although there are still some Lu-Take points lying in group A in the left bank and in group B in the right bank. As these two groups in line 1 still hint of some large amount of Lu- permeability and (group A) and Take- values (group B) the second line 2 of grouting in the targeted locations has been implemented. And, indeed, the Lu-Take values of group D for line 2 in the corresponding panel of Figure 2.6 are now coalesced closer to the origin of the Ewert’s plot, which means that in these locations the grouting curtain serves well its purpose. Moreover, its efficiency is highest in dam areas most prone to water leakage, which is where both Lu- and Take values are high, i.e. in Ewert’s group B. In fact, the three rows of panels of Figure 2.6 indicate that, when going from the exploratory check holes over line 1 to line 2 check holes, Ewerts’s group B becomes increasingly void of points, although there are still about 15% values in group A which, as stated, is difficult to grout. 2.5 Conclusions The grout curtain of Karun 4 dam is constructed on the dam foundation and abutments to prevent, or minimize the water leakage from the dam reservoir after impounding, which is the subject of the present study. Therefore, only data from the study and construction phases are used to analyse the permeability of the rock mass, before and after grouting. The measurements include Lugeon- and grout- Take values of exploratory-, line 1 and line 2 check holes of the grout curtain. The results obtained from the exploratory logs indicate the need of water proofing, whereas in special areas more treatment is necessary. Water pressure tests before and after grouting and subsequent Take values of line 1 check holes show the positions that should be grouted, which turn out to be located particularly at the entrance of the grouting galleries. Thus, here an additional second line of grouting (line 2) has been implemented. Totally, data from 64 selected boreholes are used to study the relation between Lugeon (Lu) permeability and groutability (Take) and, based on the values of these two variables, classified into four groups, following Ewert’s methodology, to infer the possible water leakage from the dam reservoir. The results indicate that the implemented grout curtains can seal leaking, as Take is lower than 50 kg/m in more than 90% of the boreholes – the requirement for efficient 30 grout Take -, although there are still areas with medium Lugeon values, i.e. data in Evert’s class A. This means that fine fissures allow water to leak, but not grout to pass. Moreover, the Lu- values indicate that the permeabilities of the right bank are generally lower than those of the left bank. This can be explained by the geological characteristics underneath the dam of the Asmari formation (AS1) subunits which show a strong correlation of Lugeon with geological zoning and existing faults in areas within the two banks, where medium permeability values are still encountered. For example, the AS1-g subunit, as key bed porous limestone has an observable effect on the water leakage and which, even after two series of grouting, leaves Lu- Take points clustered in Evert’s group A. Nevertheless, for most areas the efficiency of the two grouting curtains’ line 1 and line 2 is significant, as most of the Lu-Take points in the original exploratory hole group in the highly risky Evert’s group B (with both high Lu- and Take values) disappear in Evert’s plots for the line 1 and line 2 grout curtains. Nevertheless, there is still some scattering of Lu-Take pairs into Evert’s A- group, which means that water leakage exists at these locations. The technical solution would be to allow the extra water to seep through the drainage curtain which runs at a distance of 40 meters of the dam (see Section 2.2). The drainage curtain collects the water draining through the grout curtain and the area behind it, before being discharged to the downstream valley, and so preventing an increase of the uplift pressure which lead to rock fracture of the dam area. The seepage analysis through the drainage curtain is the topic of the ongoing study of the authors. Furthermore, as more water leakage in the bottom galleries of the left and right banks than in the two associated top galleries is to be expected after impounding, additional measurements of the piezometric heads and the effective discharges are necessary. This data will be fed in the subsequent, ongoing study, into a numerical groundwater flow model, to get a more comprehensive and quantitative picture of the water flow seepage through the dam foundation and abutments. Overall, it is to be anticipated that, as a result of all these investigations, a better evaluation of the safety and efficiency of the dam will be obtained. References Azimian A & Ajalloeian R (2013). Comparison between Lugeon with Secondary Permeability Index obtained of Water Pressure Test in Rock Masses. Electronic Journal of Geotechnical Engineering, Volume 18. 31 Barton N, Bandis S & Bakhtar K (1985). Strength, deformation and conductivity coupling of rock joints. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, Volume 22, Issue 3, pp 121–140. Berhane G & Walraevens K (2013). Geological challenges in constructing the proposed Geba dam site, northern Ethiopia. Bulletin of Engineering Geology and the Environment, Volume 72, Issue 3-4, pp 339–352. Bonacci O & Roje-Bonacci T (2008). Water losses from the Ričice reservoir built in the Dinaric karst. Engineering Geology, Voume 99, Issue 3, pp 121–127. Deere DU (1982). Cement-bentonite grouting for dams. ASCE, Geotechnical Engineering Specialty Conference on Grouting, New Orleans, pp 279–300. Ewert FK (1985). Rock Grouting with Emphasis on Dam Sites. Springer Berlin Heidelberg. Fell R, MacGregor P, Stapledon D, Bell G & Foster M (2015). Geotechnical engineering of dams. 2nd edition, CRC Press. Ford D & Williams PW (2007). Karst hydrogeology and geomorphology. [Rev. ed.]. Chichester, England, A Hoboken, NJ: John Wiley & Sons. Foyo A & Cerda J (1990). Critic permeability. New criteria for the measurement of permeability on large dams foundations. 6th International Congress International Association of Engineering Geology, Volume 2, Balkema, Netherlands, pp 1177 – 1184. Foyo A, Tomillo C, Maycotte JI & Willis P (1997). Geological features, permeability and groutability characteristics of the Zimapan Dam foundation, Hidalgo State, Mexico. Engineering Geology, Volume 46, Issue 2, pp 157–174. Foyo A, Sánchez MA & Tomillo C (2005). A proposal for a Secondary Permeability Index obtained from water pressure tests in dam foundations. Engineering Geology, volume 77, Issue 1-2, pp 69– 82. Ghafoori M, Lashkaripour GR & Tarigh Azali S (2011). Investigation of the Geological and Geotechnical Characteristics of Daroongar Dam, Northeast Iran. Geotechnical and Geological Engineering, Volume 29, Issue 6, pp 961–975. Goodman, RE, Moye DG, Van Schalkwyk A & Javandel I (1964). Ground water inflows during tunnel driving. College of Engineering, University of California. Gurocak Z & Alemdag S (2012). Assessment of permeability and injection depth at the Atasu dam site (Turkey) based on experimental and numerical analyses. Bulletin of Engineering Geology and the Environment, Volume 71, Issue 2, pp 221–229. Haghnejad A, Ahangari K & Noorzad A (2014). Investigation on Various Relations Between Uniaxial Compressive Strength, Elasticity and Deformation Modulus of Asmari Formation in Iran. Arabian Journal for Science and Engineering, Volume 39, Issue 4, pp 2677–2682. Hosseiny Sohi SM, Koch M & Ashjari J (2014). Construction Effects of the Karun 4 Dam, Iran, on the Groundwater in the Adjacent Karstic Aquifer. ICHE 2014, Hamburg - Lehfeldt & Kopmann (eds). 32 Houlsby AC (1990). Construction and design of cement grouting. A guide to grouting in rock foundations. New York, Wiley (Wiley series of practical construction guides). Jaeger Ch (1979). Rock mechanics and engineering. Cambridge University Press. Kocbay A & Kilic R (2006). Engineering geological assessment of the Obruk dam site (Corum, Turkey). Engineering Geology, Volume 87, Issue 3, pp 141–148. Kutzner Ch (1996). Grouting of soil and rock. Rotterdam: Balkema. Lugeon M (1933). Barrages et Geologie. Dunod, Paris Mahab Ghods Consulting Engineering Co. (2002). Final engineering geological report (phase 2). Karun 4 dam project. Mahab Ghods Consulting Engineering Co. (2010). Final excavation and injection report of the Karun 4 dam grout curtain. Milanovic P (2004). Water resources engineering in Karst. CRC Press. Nonveiller E (1989). Grouting, theory and practice. Amsterdam, New York: Elsevier (Developments in geotechnical engineering, 57). Rajabi AM Khodaparast M & Edalat A (2015). Investigation of the geological and geotechnical characteristics of the Tanguyeh dam site in southeastern Iran. Bulletin of Engineering Geology and the Environment, Volume 74, Issue 3, pp 861–872. Sadeghiyeh SM, Hashemi M & Ajalloeian R (2013). Comparison of Permeability and Groutability of Ostur Dam Site Rock Mass for Grout Curtain Design. Rock mechanics and rock engineering, volume 46, Issue 2, pp 341–357. Singh B & Goel RK (1999). Rock mass classification. A practical approach in civil engineering. 1st ed. Amsterdam, New York: Elsevier. Turkmen S, Özgüler E, Taga H & Karaogullarindan Talip (2002). Seepage problems in the karstic limestone foundation of the Kalecik Dam (south Turkey). Engineering Geology, Volume 63, Issue 3, pp 247–257. Uromeihy A (2000). The Lar Dam; an example of infrastructural development in a geologically active karstic region. Journal of Asian Earth Sciences, Volume 18, Issue 1, pp 25–31. Uromeihy A & Barzegari G (2007). Evaluation and treatment of seepage problems at Chapar-Abad Dam, Iran. Engineering Geology, Volume 91, Issue 2, pp 219–228. Uromeihy A & Farrokhi R (2012). Evaluating groutability at the Kamal-Saleh Dam based on Lugeon tests. Bullten of Engineering Geology Environmental, Volume 71, Issue 2, pp 215–219. Yang, ChP (2004). Estimating cement take and grout efficiency on foundation improvement for Li-Yu- Tan dam. Engineering Geology, Volume 75, Issue 6, pp 1–14. 33 Chapter 3- Monitoring of Seepage Controlling System of Karun 4 Dam in Iran Abstract Karun 4 concrete dam with 230 m height is the first dam constructed at the upstream of Karun river in south west of Iran. The dam site area is located on the carbonate limestone of the Asmari formation which has strong karstic characteristics and, subsequently, a high potential for water leakage. To deal with this problem, a seepage controlling system is implemented within the dam that includes a grout curtain as a virtual barrier against water leakage on the dam foundation and abutments, and a drainage curtain at the downstream of the dam body to allow collection of the escaped water from the grout curtain; and the monitoring instruments: Standpipe piezometers and observation wells. The object of present study is to present the dam seepage controlling system; and to investigate the water head variation of the piezometers and observation wells, using obtained measurements after impounding, to evaluate the efficiency of the grout curtain. In addition, the areas which have larger water head variations are discussed, indicating the variations are generally higher in the left- than in the right bank. And finally, the need of a numerical model of the groundwater flow of the whole area to project the dam reservoir behavior for the coming years, is concluded. Key words: Seepage, Monitoring, Karun 4 dam, Stand-pipe piezometer, Observation well, grout curtain, drainage curtain, carbonate aquifer, Karst. 34 3.1 Introduction It is generally agreed upon today that the monitoring of a large dam plays a big role in the operation phase, particularly, in its first years. The dam’s stability and safety can only be assured through the continuous observation of different parameters, such as displacements, pore-pressure, temperature, joints and cracks, water seepage, followed by documentation and thorough analysis of the obtained data. Water seepage through the dam abutments and foundation can be prevented or reduced via the construction of a grout curtain (Houlsby,1990; Fell et al. 2015). Then the excess water can be discharged downstream through a drainage curtain. This whole process can be monitored by installed piezometers and observation wells (Bartholomew & Haverland, 1978). Seepage controlling system are of particular importance for dams constructed within porous limestone formations, because in such areas an increase of the water level in the reservoir may lead to a partial dissolution of the carbonate formation, a process known as karstification (Milanović, 2004). The Lar dam in Iran suffers from this problem and for this reason its reservoir was never fully impounded (Uromiyehi 2000). Turkmen et al. (2002) evaluated the seepage problems in the karstic limestone foundations of the Kalecil dam in Turkey, whereas Ghobadi et al. (2005) studied the seepage through the karstic limestone of the right abutments of the Karun 1 dam in Iran. Kocbay & Kilic (2006) assessed the Obruk dam site in Turkey from a geological engineering point of view. Mohammadi et al. (2007) presented a method for leakage study in the karst foundations of the Khersan 3 dam, Iran, while Fazeli (2007) studied the construction of a grout curtain in the karst environment of the Salman Farsi Dam, Iran. Uromeihy & Barzegari (2007) studied the seepage problem of the Chapar-Abad dam in Iran with water pressure tests and seepage modelling, using the PLAXIS software. The water losses of the Ricice reservoir built in the Dinaric karst in Croatia have been studied by Bonacci and Roje-Bonacci (2008). Chen et al. (2010) presented a new classification of seepage controlling mechanisms in geotechnical engineering. Ghafoori et al. (2011) studied the Daroongar dam in Iran from a geological and geotechnical point of view, whereas Dafny et al. (2015) evaluated the temporal changes in hydraulic conductivities in the karst terrain around the Dokan dam in Kurdistan Iraq. 35 In the present paper, the seepage controlling system of the Karun 4 dam in Iran is described and studied. This dam is constructed on the Asmari formation which includes porous limestone with a high potential of karstification in some of its subunits. Therefore, an expanded grout curtain in the form of different grouting galleries has been constructed in the dam abutments and foundation. Water head changes downstream of the reservoir are monitored based on measurements in piezometers installed in the galleries of the grout curtain and in observation wells placed at different points of the dam site from the grouting galleries to the access- and drainage galleries downstream of the dam. These measured water heads are used to study the groundwater situation around the dam area after impounding and in the operation phase. 3.2 Seepage controlling system of the Karun 4 dam 3.2.1 Site specification and geology The Karun 4 dam is located close to the source of the Karun River in the state of Chaharmahal and Bakhtiari, in southwest Iran (Figure 3.1). Its purpose is to regulate the Karun river’s discharge interconnected with other under-operation dams, such as Karun 3-, Karun 1-, Masjed-Soleiman- and Gotvand dams, and to control the frequent destructive floods of the Karun River, but also to generate about 1GW hydroelectric energy (Table 3.1). The dam has been impounded since March 2010. Table 3.1. Karun 4 dam technical specifications Karun 4 concrete double arc dam specifications Specifications Height from foundation 230m Crest length / width 440m / 7m Total volume of the reservoir 2190 million m3 Area of reservoir 29 km2 Normal water level 1025 m.a.s.l. Turbines capacity of powerhouse 4 x 250 MW Hydroelectric g eneration 2100 GWh/year The Karun 4 dam is located in the tectonically folded belt of the Zagros Mountains within a canyon located at the southwestern flank of the Kuh-Sefid anticline which has a NW-SE trend. Most of the reservoir is situated on Cretaceous to Miocene limestone and marly limestone (Haghnejad et al., 2014). 36 Figure 3.1. Location of Karun 4 dam in Iran with observation wells’ layout The geological formations around the dam site consist of carbonate layers of the Asmari formation (As) which can be divided into lower (As1), middle (As2) and upper (As3) units and marlstone and marly limestone of the impervious Pabdeh formation (Pd) upstream of the dam site (see Figure 3.1). The geology of the dam site is affected by 11 faults in the left abutment and 9 faults in right abutment, all of which are mainly dextral with their dips in the range of 25°-85°. The joint openings along the bedding at the ground surface range from 1 to 100 mm and decrease with increasing depth. The lower unit of the Asmari formation AS1, consisting of limestone, marly limestone and marlstone with some porous limestone and with an apparent thickness of about 200 m in normal direction to the bedding, constitutes the abutments and foundation rock of the dam (Hosseiny Sohi et al., 2014). These units outcrop in the dam valley abutments and are highlighted in the cliff around the site. The strikes of the beds are nearly parallel to the dam axis and dip at angles of 70-80 degrees towards the downstream. The AS1 unit is subdivided into 10 subunits, named AS1-a to AS1-j, from upstream to downstream. The bedding thickness of the rock mass varies from 2 m to 30 m. A porous limestone layer (AS1-g) was 37 Figure 3.2. Grout curtain expanded section view and the piezometers’ layout identified, which had a vuggy feature and was recognized as a key bed underneath the dam site. The karstification degree decreases gradually below the river’s water level. 3.2.2 Grout curtain To seal the dam foundation and abutments, a very large grout curtain was designed and constructed. The construction was done through the excavation of 5 serial galleries in both left and right abutments at different elevations from the dam crest elevation to the bottom, ending with the foundation gallery # 806 (see Figure 3.2). From inside these galleries, the grouting boreholes are drilled, examined and injected, to have an integrated virtual wall or barrier, called grout curtain against water leakage. The drill of the grouting holes was done using the split- spacing boreholes method (Houlsby, 1990), wherefore from the starting point of each gallery (kilometer 0+00, the intersection of the concrete dam body and excavated gallery inside the rock mass), every 16 meter the primary holes or A series were drilled, the secondary holes, B series, were done between them at intervals of 8 meters, the tertiary holes, C series, every 4 meters, the quaternary holes, D series, every 2 meters and, finally, the quinary holes, E series, every meter. In fact, as the primary holes were the first drilled boreholes, they are used as a reliable source of the original in situ data representing the existing mechanical properties of the rock mass. Some of grouting holes were selected as exploratory boreholes and rock permeability (Luegon) tests and geological logs were carried out in them. They were then considered in the grouting phase as a normal drilled borehole and the amounts of their cement Take were reported. Details and results of these analyses can be found in Hosseiny Sohi et al. (2017). Based on the obtained 38 results from these exploratory boreholes, the grout curtain in two (in some special areas in three) lines was designed and constructed, in order to sew/stitch the dam structure to the Pabdeh impermeable formation upstream. The first line was done almost along the galleries, whereas the second line located 0.9 m upstream toward the reservoir of the former, was drilled principally in the initial sections of the galleries, because, after grouting, the farther sections the first line of galleries exhibited very few permeability and there was thus no need to grout them in addition. Finally, after performing the second line, the third line was drilled and grouted only in those areas in the lower galleries which still showed more permeability than the allowable limit. After implementation of the first and second lines of the grout curtain, the check holes were drilled and the permeability and Take were measured (see Hosseiny Sohi et al., 2017, for further details). 3.2.3 Drainage curtain Due to the mentioned geological characteristics of the Asmari formation rock and based on the measured data from the check holes, it is evident that the grout curtain constructed is not able to completely seal off the foundation and abutments of the dam, so that remaining seepage water must be properly drained. Therefore, an additional drainage curtain was implemented whose purpose is not only to detect and to explore the still leaking water through the grout curtain, but also to reduce the uplift pressure due to piping and so to prevent the breaking of the rock mass downstream. This grout curtain consists of boreholes drilled at a horizontal distance of 40 meters in downstream direction in two galleries in the left (LDG1 and LDG2) and right abutment (RDG1 and RDG2) (see Figure 3.3). The elevation of LDG1 and RDG1 is +932 m.a.s.l, whereas the LDG2 elevation increases from 853 m at the beginning to 867.27 m at the end, and RDG2 from 860 m to 875.11 m. In the first 70 meter of the beginning of all four galleries, holes are drilled every 6 meter, and every 12 meter afterwards. The depths of the top galleries’ holes (LDG1 and RDG1) range from 50m to 71m, whereas in the bottom galleries (LDG2 and RDG2) the former are 100m at the first 70m of the gallery and 50m afterwards (Mahab Ghods report, 2010). The amount of drained water is measured in the determined points, especially, at the final sump of the bottom drainage galleries (LDG2 and RDG2), where the water is discharged naturally to the downstream valley. 39 Figure 3.3. Cross-sectional view of the grout curtain system 3.2.4 Monitoring instruments To monitor the dam’s geotechnical and hydraulic behaviour, a network of instruments is installed in the dam body, grouting galleries and the different stations of the dam site. The most important instruments of the seepage controlling system are the standpipe piezometers and observation wells. The former, also known as Casagrande piezometers, are installed in the grouting galleries upstream and downstream of the grout curtain (see Figure 3.2), whereas the observation wells are drilled in the different points of the dam body, abutments, the drainage curtain and downstream of the dam site (see Figure 3.1). 3.2.4.1 Standpipe piezometers The Casagrande standpipe piezometers are used to measure the hydraulics/piezometric heads (water levels) in boreholes, to monitor the effectiveness of the grout curtain and to observe the leakage and the groundwater movements in the rock abutments of the dam. They are comprised of two parts: at its lowest point is a porous piezometer tip; connected to the tip is a riser pipe which continues upwards out of the top of the borehole. To measure the water level, the filter tip zone is packed with sand and then backfilled above. Three different types depending on the location of the piezometer sensor are distinguished: Type 1-upstream of grout curtain, type 2- between grout and drainage curtain and type 3- downstream of drainage curtain. According to 40 the monitoring report of Mahab Ghods engineering company (2015), the heads in the standpipe piezometers at the dam site are measured daily. 3.2.4.2 Observation wells The observation wells are drilled in the grouting galleries, drainage galleries and access gallery in order to monitor the groundwater levels of the dam abutments and foundation after impounding the reservoir and during the operation phase. The unlined observation wells, unlike the lined standpipe piezometers, return the free-standing groundwater table. The wells drilled in the grouting galleries have an inclination (from horizontal) of 70 degrees, whereas the wells downstream of the grout curtain are drilled vertically (the observation well at the first left drainage gallery is an exception with a 70-degree inclination). 3.3 Data Mahab Ghods consulting engineering company provided geology maps, layout, view and details of the grout curtain maps, the engineering geology report “Final excavation and injection report of the Karun 4 dam grout curtain”, the rock mechanics report, and the geotechnical and geological data of the drilled holes such as geological logs, and water pressure test results (Lugeon values and grout take), wherefore the latter have been analysed and evaluated by Hosseiny Sohi et al. (2017). In the present study, data from the first impounding in March 2010 until November 2014 are used, during which time period the standpipe piezometers were read out daily in the grouting galleries of the left and right banks, whereas the observation wells were sampled on a monthly interval. In each bank, 16 piezometers at the downstream side are considered, however, since in the right bank all four piezometers of the RG2 and the two first piezometers of the RG3 are dry, the results of a total of 26 piezometers are used. For the groundwater levels, data of 11 observation wells are analysed. With this data, together with the grout curtain results before impounding analysed by Hosseiny Sohi et al. (2017), it is possible to obtain a comprehensive picture of the seepage controlling system of the dam. 3.4 Results 3.4.1 Standpipe piezometers The data of the standpipe piezometers which are installed in all grouting galleries of the left and right banks, except the first right gallery (RG1), originate only from the piezometers 41 located downstream of the grout curtain (type 2, see Section 3.2.4.1). This is because, the piezometers of the first left gallery (LG1), the second right gallery (RG2) and two first piezometers at the third right gallery (RG3) downstream of the grout curtain were during the measuring period mentioned. The time-series graphs of the piezometric heads of the third, fourth and fifth lower rows of the grout curtain galleries are illustrated in Figure 3.4. In each grouting gallery, four piezometers at the different distances from the dam body have been sampled with the results summarized in Table 3.2. Furthermore, the water head variations for four piezometers for each gallery are shown for the left and right bank in the two bar plots of Figure 3.5, respectively. Figure 3.5 shows that, despite a few anomalies in the piezometric time series, there is no critical trend of an increasing water head; for instance, in LG3, all water levels are fluctuating within a narrow range. Nevertheless, at the same elevation at the right bank, the head differences between the only two active piezometers are big. Moreover, the variations of the piezometric heads for most cases are essentially a function of the reservoir’s water level at the top, which is indicated at the right axis of the time series panels. Thus, the piezometer levels are increasing in the first months of impounding, after which time they reach their corresponding, almost constant elevation, which is function of the seasonal changes and of the reservoir water level. An exception is the fourth piezometer SP250 of LG5 which shows irregular results with the highest head variations. Moreover, whereas in the left bank the water heads increase when going from galleries at higher elevations to those at lower elevations, this trend is not observable at the right bank where the highest values belong to RG4 and not to RG5. Also, except for RG4 and RG5, in all other galleries, the water head variations increase along the gallery from the first piezometer toward the fourth one. At RG4 the second piezometer SP100 has lower values than the first one, SP50, while for RG5 the second piezometer SP100 has the highest head value in that group. However, the overall maximally measured head values belong to SP280 of RG4. 42 Figure 3.4. Water heads measured by the downstream piezometers on left and right banks of the grout curtain of the 3rd gallery (top row), 4th gallery (middle row) and 5th gallery (bottom row) 840 860 880 900 920 940 960 980 1000 1020 1040 900 920 940 960 980 1000 March-10 March-11 March-12 March-13 March-14 March-15 P ie z o m e tr ic E le v a ti o n ( m ) Piezometers of Gallery LG3 - Left Bank SP 50 SP 100 SP 170 SP 233 R.W.L. R e s e rv o ir E le v a ti o n ( m ) Date 840 860 880 900 920 940 960 980 1000 1020 1040 900 920 940 960 980 1000 March-10 March-11 March-12 March-13 March-14 March-15 P ie z o m e tr ic E le v a ti o n ( m ) Piezometers of Gallery RG3- Right Bank SP 185 SP 265 R.W.L. R e s e rv o ir E le v a ti o n ( m ) Date 840 860 880 900 920 940 960 980 1000 1020 1040 860 880 900 920 940 960 980 1000 March-10 March-11 March-12 March-13 March-14 March-15 P ie z o m e tr ic E le v a ti o n ( m ) Piezometers of Gallery LG4 - Left Bank SP 50 SP 105 SP 148 SP 250 R.W.L. R e s e rv o ir E le v a ti o n ( m ) Date 840 860 880 900 920 940 960 980 1000 1020 1040 860 880 900 920 940 960 980 1000 March-10 March-11 March-12 March-13 March-14 March-15 P ie z o m e tr ic E le v a ti o n ( m ) Piezometers of Gallery RG4 - Right Bank SP 50 SP 100 SP 200 SP 280 R.W.L. R e s e rv o ir E le v a ti o n ( m ) Date 840 860 880 900 920 940 960 980 1000 1020 1040 840 860 880 900 920 940 960 March-10 March-11 March-12 March-13 March-14 March-15 P ie z o m e tr ic E le v a ti o n ( m ) Piezometers of Gallery LG5 - Left Bank SP 50 SP 100 SP 170 SP 250 R.W.L. R e s e rv o ir E le v a ti o n ( m ) Date 840 860 880 900 920 940 960 980 1000 1020 1040 840 860 880 900 920 940 960 March-10 March-11 March-12 March-13 March-14 March-15 P ie z o m e tr ic E le v a ti o n ( m ) Piezometers of Gallery RG5 - Right Bank SP 50 SP 100 SP 170 SP 250 R.W.L. R e s e rv o ir E le v a ti o n ( m ) Date 43 Figure 3.5. Water head variations of the four piezometers of each gallery in left and right bank Table 3.2. Standpipe piezometers results of the left and right banks of the grout curtain In summary, it can be deduced from the visual inspection of Figure 3.5 that the head variations in the left bank are larger than those in the right bank. This would confirm the results of a previous study of the authors (Hosseiny Sohi et al. 2017) who determined the (Lugeon) permeability of the grout curtain of the Karun 4 dam. On the other hand, no reasonable relation between the observed water head variations and the geology of the dam site can be deduced; thus the water pressure tests results obtained from the exploratory boreholes of the grouting curtain (series A boreholes) (see Hosseiny Sohi et al. 2017) return no significant trend that could be related to the piezometers results; for example the Lugeon value at the sensor elevation of the fourth piezometer SP250 of the LG5 is 23 Lu, whereas in many other positions the Lugeon values go up to more than 60 Lu. On the other 0 10 20 30 40 50 60 70 80 90 LG2 LG3 LG4 LG5 W a te r H e a d V a ri a ti o n ( m ) Water Head Variation of Piezometers (Left Bank) 1st Piezometer 2nd Piezometer 3rd Piezometer 4th Piezometer 50 100 170 233 LG3 50 105 148 250 LG4 50 100 170 250 LG5 40 100 180 250 LG2 0 10 20 30 40 50 60 70 80 90 RG3 RG4 RG5 W a te r H e a d V a ri a ti o n ( m ) Water Head Variation of Piezometers (Right Bank) 1st Piezometer 2nd Piezometer 3rd Piezometer 4th Piezometer 50 115 185 265 RG3 50 115 200 280 RG4 50 100 170 250 RG5 Left Galleries Sensor Elevation Distance from dam(m) Mean (m) Water Head Variation (m) Right Galleries Sensor Elevation Distance from dam(m) Mean (m) Water Head Variation (m) 40 952,7 5,8 36 Dry - LG2 100 952,4 4,7 RG2 65 Dry - (974 m) 180 952,1 7,0 (997 m) 130 Dry - 250 950,7 5,5 210 Dry - 50 910,9 2,6 50 Dry - LG3 100 912,3 5,6 RG3 115 Dry - (932 m) 170 914,7 14,6 (932 m) 185 929,7 28,6 233 915,9 11,2 265 957,7 23,5 50 879,0 27,4 50 871,3 7,7 LG4 105.5 883,9 26,9 RG4 115 868,3 4,9 (890 m) 148.5 891,1 28,6 (890 m) 200 879,8 21,6 250 916,2 46,7 280 945,7 83,1 50 876,7 31,6 50 870,8 42,8 LG5 100 881,7 45,9 RG5 100 881,4 53,3 (848 m) 170 903,7 47,9 (848 m) 170 869,6 42,2 250 910,1 74,5 250 869,7 39,1 810m 865m 965 m 910 m 950 m 910 m 865m 810m 44 hand, for the corresponding fourth piezometer in LG4 the Lugeon value is only 11 Lu.The reasons for these irregular changes may be due to the grout curtain and, probably, it’s construction, since after carrying out the operation of the cement injection through the grouting galleries, there is still water leakage at the dam site. To evaluate entirely this process and having a general view of the water flow from the reservoir downstream through the rock mass abutments of the dam, a comprehensive numerical model is needed, which is the issue of the ongoing study of the authors. 3.4.2 Observation wells Table 3.3 lists the groundwater elevations measure in the observation wells which are drilled in the grouting galleries of both banks and also in the drainage galleries and the corresponding access galleries of the dam site. The observation wells names indicate the observations well of the left bank (LOW1 to LOW8) and of the right bank (ROW1 to ROW3). Three wells in LG3 and two wells in LG4 are located in the grouting galleries, whereas all other six wells are at different distances from the dam abutments and grout curtain. The drilling dip of the wells in the grouting galleries and of well LOW6 in the first left drainage gallery is 70 degrees from the horizontal in downstream direction; meanwhile all other wells far from the grout curtain are drilled vertically. At the left bank, LOW1 is the farthest from the grout curtain (225 meter) and at the right bank, ROW2 is 200 meters away from the grout curtain. Table 3.3. Results of observation wells’ measurements ----------------------------------------------------------------------------------------------------------- Top Bottom Dip Elevation (m) Elevation (m) (Degree) Mean (m) Variation (m) LG2 LOW1 963 840 Vertical 872.0 18.8 LOW2 932 888 70 915.4 18.8 LOW3 932 871 70 883.3 5.5 LOW4 932 871 70 886.5 9.1 LOW5 932 832 Vertical 864.4 27.9 LOW6 932 838 70 863.6 12.1 LOW7 893 857 70 882.7 5.2 LOW8 893 841 70 883.3 5.2 ROW1 927 858 Vertical 871.6 23.6 ROW2 931 847 Vertical 865.3 18.5 RG4 ROW3 884 840 Vertical 849.8 12.9 LG3 LG4 RG3 GroundwaterObservation Well Gallery 45 Figure 3.6. Groundwater level measurements of the observation wells The time series of all observation well measured water elevations are illustrated in Figure 3.6. Since there are only three wells in the right bank, the data of both banks are drawn in one figure and are not separated. Figure 3.6 and, more quantitatively, Table 3.3, indicates that the largest temporal changes of the groundwater level happen in the left bank well LOW5 which has also the deepest bottom elevation of 832 meter. This well is located at the access gallery of LG3, at a distance of 85 meters from the grout curtain. At this bank the lowest variation belongs to the LOW7 and LOW8, both at the grouting gallery LG4. At the right bank ROW1 at a distance of 82 meters from RG3 shows the highest ground water elevation variations. Generally, the variations of the ground water levels are larger in the left than in the right bank. This finding agrees with that obtained from the standpipe piezometers results in Section 3.4.1 and confirms the achievements of the grout curtain study of the dam before impounding, done by Hosseiny Sohi et al. (2017). It is also in compliance with the geological characteristics of the Asmari formation on the left bank and the water pressure test results after two series of cement injection at the galleries of the grout curtain. 820 840 860 880 900 920 940 960 980 1000 1020 820 840 860 880 900 920 940 960 March-10 March-11 March-12 March-13 March-14 March-15 G ro u n d w a te r L e v e l ( m ) LOW 1 LOW 2 LOW 3 LOW 4 LOW 5 LOW 6 LOW 7 LOW 8 ROW 1 ROW 2 ROW 3 R.W.L. R e s e rv o ir E le v a ti o n (m ) Date 46 3.5 Conclusions Monitoring of large dams is one of the most important tasks for ensuring the dam’s stability and proper maintenance after construction of the dam, especially, in the first years of the operation phase. The Karun 4 dam, as the highest dam of Iran and as a young dam, needs particular monitoring attention, namely, of seepage and stability. Because of the presence of the Asmari limestone formation underneath the dam foundation and abutments, the karstification is of most concern to this project, as this could lead to water leakage from the reservoir. Therefore, a comprehensive seepage controlling system with the following features in is established: Through the grouting galleries, an expanded grout curtain is drilled and injected with cement grout in the dam foundation and both abutments, to prevent the water leakage. Since the grout curtain is not able to absolutely seal the water passage, a drainage curtain is constructed at distance of 40 meter further downstream to discharge the seeping water downstream in the river valley. This drainage curtain serves also to prevent the increase of uplift pressure at the dam foundation and abutments. To observe the efficiency of the grout curtain, standpipe piezometers are installed in the grouting galleries; whereas the observation wells for the measurement of the groundwater levels are located in the same galleries and at other places further downstream. The results of the standpipe piezometers head measurement indicate that the grout curtain works generally well, however, there are some piezometers which show large head variations, which are generally higher in the left than in the right bank. Furthermore, it the left bank, the water heads increase proportionally with the depth of the piezometer to that they reach their maximum values at the bottom gallery (LG5). The groundwater level changes measured by the observation wells confirm these results, such that the changes in the left bank are larger than in the right bank. With the data available so far, is not possible to project the behaviour of the grout curtain and of the reservoir for the coming years. Therefore, a numerical model of the groundwater flow in the whole area, taking into account the geological characteristics of the dam site, and using the observed measurements of the monitoring instruments for calibration is recommended, which is the issue of the ongoing study of the authors. In comparison with other technics such as dye- tracing tests, which have already been done (6 series) at the dam site, numerical modelling has not only an economic justification, but allows also to simulate various scenarios; thus, giving 47 the possibility to predict the future trends of the groundwater movements underneath the Karun 4 dam. References Bartholomew ChL, & Haverland ML (1987). Concrete dam instrumentation manual. United Stated Department of the Interior Bureau of Reclamation (USBR). Bonacci O, & Roje-Bonacci T (2008). Water losses from the Ričice reservoir built in the Dinaric karst. Engineering Geology, Voume 99, Issue 3, pp 121–127. Chen Y, Hu R, Zhou Ch, Li D, Rong G, & Jiang Q (2010). A new classification of seepage control mechanisms in geotechnical engineering. Journal of Rock Mechanics and Geotechnical Engineering, Volume 2, issue 3, pp 209–222. Dafny E, Tawfeeq KJ, & Ghabraie K (2015). Evaluating temporal changes in hydraulic conductivities near karst-terrain dams: Dokan Dam (Kurdistan-Iraq). Journal of hydrology, Volume 529, pp 265– 275. Fazeli MA (2007). Construction of grout curtain in karstic environment case study: Salman Farsi Dam. Environmental geology, Volume 51, Issue 5, pp 791–796. Fell R, MacGregor P, Stapledon D, Bell G, & Foster M (2015). Geotechnical engineering of dams. 2nd edition. Ghafoori M, Lashkaripour GR, & Tarigh Azali S (2011). Investigation of the Geological and Geotechnical Characteristics of Daroongar Dam, Northeast Iran. Geotechnical and Geological Engineering, Volume 29, Issue 6, pp 961–975. Ghobadi MH, Khanlari GR & Djalaly H (2005). Seepage problems in the right abutment of the Shahid Abbaspour dam, southern Iran. Engineering Geology, Volume 82, Issue 2, pp 119–126. Haghnejad A, Ahangari K & Noorzad A (2014). Investigation on Various Relations Between Uniaxial Compressive Strength, Elasticity and Deformation Modulus of Asmari Formation in Iran. Arabian Journal for Science and Engineering, Volume 39, Issue 4, pp 2677–2682. Hosseiny Sohi SM, Koch M & Ashjari J (2014). Construction Effects of the Karun 4 Dam, Iran, on the Groundwater in the Adjacent Karstic Aquifer. ICHE 2014, Hamburg - Lehfeldt & Kopmann (eds) Hosseiny Sohi SM, Koch M & Ashjari J (2017). Evaluating permeability and groutability at Karun 4 dam Iran using Lugeon values and grout Take. The 85th Annual Meeting of International Comission on Large Dams, Prague, July 3-7, 2017. Houlsby AC (1990). Construction and design of cement grouting. A guide to grouting in rock foundations. New York, Wiley (Wiley series of practical construction guides). Kocbay A & Kilic R (2006). Engineering geological assessment of the Obruk dam site (Corum, Turkey). Engineering Geology, Volume 87, Issue 3, pp 141–148. 48 Mahab Ghods Consulting Engineering Co. (2002). Final engineering geological report (phase 2). Karun 4 dam project. Mahab Ghods Consulting Engineering Co. (2010). Final excavation and injection report of the Karun 4 dam grout curtain. Mahab Ghods Consulting Engineering Co. (2015). Monitoring report of Karun 4 dam. Milanovic P (2004). Water resources engineering in Karst. CRC Press. Mohammadi Z, Raeisi E & Bakalowicz M (2007). Method of leakage study at the karst dam site. A case study: Khersan 3 Dam, Iran. Environmental geology, Volume 52, Issue 6, pp 1053–1065. Turkmen S, Özgüler E, Taga H & Karaogullarindan T (2002). Seepage problems in the karstic limestone foundation of the Kalecik Dam (south Turkey). Engineering Geology, Volume 63, Issue 3, pp 247– 257. Uromeihy A (2000). The Lar Dam; an example of infrastructural development in a geologically active karstic region. Journal of Asian Earth Sciences, Volume 18, Issue 1, pp 25–31. Uromeihy A & Barzegari G (2007). Evaluation and treatment of seepage problems at Chapar-Abad Dam, Iran. Engineering Geology, Volume 91, Issue 2, pp 219–228. 49 Chapter 4 - Simulation of Groundwater Flow in the Dual Porous Media of the Karun 4 Dam Foundation and Abutments using Conduit Flow Process (CFP) for MODFLOW-2005 Abstract Water leakage through dam foundation and abutments is one of the prevalent problems when a dam is constructed in a carbonate aquifer with karst potential. The Karun 4 dam in Iran is built in a porous limestone formation that could lead to water leakage. The geological field investigations before its construction indicated the zones with fractures and conduits. To mitigate this problem, a large grout curtain has been drilled and injected with cement grout via different grouting galleries in the dam foundation and abutments. To evaluate the function of the curtain, and find out the possible water leakage paths, a numerical model of dual porous media is set up using Conduit Flow Process (CFP) for MODFLOW-2005. The model simulates both porous matrix laminar- and conduits laminar or turbulent flow. Using the reservoir water levels, drainage observations, observed groundwater head from installed piezometers and observation wells, the model is successfully calibrated and validated. In the transient state, firstly, the matrix medium is calibrated in MODFLOW as a laminar groundwater flow; then, the conduits are set in the CFP model. The conduit layout is determined by examination of the faults from geological maps. The approved faults are assigned as pipes in the model. The pipes calibration parameters are their hydraulic conductivity and diameter. Afterwards the efficiency of the grout curtain is evaluated which identifies various hydraulic conductivities along the grout curtain. The present model enables one to investigate different scenarios of water leakage in the model area which is the subject of a subsequent study. Key words: Groundwater flow, Carbonate aquifer, Karst, Dual porous media, Numerical modelling, CFP, MODFLOW, Water leakage, Karun 4 dam, Iran 50 4.1 Introduction During dam construction and later in the operating phase, various structural, geotechnical and hydrogeological problems can take place. One of the most challenging hydrogeological problems is water leakage. Leakage not only causes water loss of the reservoir, but can, at the same time, threaten the dam stability. A well-known example of excessive reservoir leakage over the years is the Montejaque dam in Spain. Montejaque dam has yet to be filled to its maximum capacity with water after many years of construction, because of reservoir and abutments leakages (Sahuquillo,1985). Another example is the Camarassa dam in Spain, where the original 2.5 m³/s leakage has increased to 11.5 m³/s over the years (Jaeger,1979). The May dam in Turkey that was constructed on karst limestone has never been filled (Doğan and Çiçek, 2002). A final example is the Lar dam in Iran which has 10.5 m³/s leakage nowadays (Uromeihy,2000). One of the potentially paths for water leakage in dams is fractured and faulted zones, which could conduct water flow from upstream to downstream (Knill ,1972). Flow theory in a classical porous media was presented by Darcy (1856) based on experimental studies he did. But formulating fluid movement theory in fractured and faulted zones has been facing a lot of problems, because of the different sizes of empty spaces between particles which precludes the use of classical porous media flow theory. Classification of karst aquifers from the porosity point of view are: 1- Primary porosity which is established during sedimentation and diagenesis. Karst- aquifers store only small amounts of water, normally. 2- Secondary porosity describes fractured and faulted zones which have been formed by tectonics forces. 3- Third porosity is that which has been formed by solution processes (Worthington et al., 2000). Ford and Williams (2007) discussed the multiplicity of porosity that causes complicated flow regimes in karsts in more detail. A lack of precious data of conduit geometry has intensified such complications (Mohrlock and Teutsch, 1977). Quinlan and Ewers (1985) have classified karst aquifers based on their hydraulic specifications: 51 1- Mature system with a conduit system that has a high hydraulic conductivity (e.g., Alps karst system), which has been formed by fracturing and faulting processes or by the dissolution of conduits. In such a media groundwater flow is essentially a turbulent pipe flow, as the pore openings are more than 1 cm wide. 2- Dispersion system which has a conduit matrix with openings less than 1cm wide. Groundwater tends to behave as laminar Darcy flow. 3- Hybrid system which is a combination of the above two mentioned systems. In this case both conduit and matrix systems co-exist. Most karst aquifers can be classified as a hybrid system, whereby the karst features lead to an aquifer with high flow velocities (high hydraulic conductivities) and significant storage of water. Therefore, any attempt of karst aquifer modeling must take into account the special porosity structure of such a media, namely, the so-called double-porosity continuum. This duality concept was first presented by Barenblatt et al. (1960). It was developed further and its applicability to Karstic aquifers proved by Atkinson (1977), Gunn (1985), Yurtsever (1986), Huntoon (1995), Zhang (1996), Keeler (1997), Kiraly (1998), Halihan & Wicks (1998), White (1999), Jeannin (2001), Neumann (2005), and Ghasemizadeh et al. (2012). More recent numerical models developed for groundwater flow simulations in a karst media are: - Double Continuum Porous Equivalent model by Sauter (1992), - CAVE (Carbonate Aquifer Void Evolution) model by Clements et al. (1996), - GW (Groundwater model) by Cornaton and Perroche (2002), and - CFP (Conduit Flow Process) for MODFLOW-2005 by Shoemaker et al. (2008). In all these numerical approaches the matrix media is considered separately from the conduit media. The two medias are connected hydraulically by a coefficient named the exchange coefficient. The most important difference between the models mentioned is how the conduit medium is defined. Thus, in the Double Continuum Porous Equivalent model, the conduit medium is considered by a hypothetical and equivalent porous medium, where the Darcy linear flow law is applicable. In the CAVE and CFP models, the conduits are presented by a pipe network, in which turbulent flow occurs. The GW model is essentially like the CAVE- and CFP models, but, instead of a finite difference method, a finite element method is used. 52 The most widely used codes in the hydrogeological community are: the US Geological Survey codes of MODFLOW, starting from MODFLOW 88 (McDonald and Harbaugh, 1988), MODFLOW-2000 (Harbaugh et al., 2000) to MODFLOW-2005 (Harbaugh, 2005). Thus, CFP is applied in works of many authors like Hill et al. (2008), Reimann et Hill (2009); Kuniansky and Shoemaker (2009), Ashok and Sophocleous (2009), Reimann et al. (2011); Saller et al. (2013), Gallegos et al. (2013), and Assari and Mohammadi (2017). Hartmann et al. (2014) provided a remarkable summary of the hydrological approaches in Karst water resources. They explored different conceptual models of karst systems and how they can be translated into numerical models;discussed the limitations of current karst models, and concluded by providing new research directions. The present study relies on the results of the recent works of the authors regarding Karun 4 dam project: - Chemical investigations of water in dam reservoir and downstream that proves the significant arising of the ground water level in the area (Hosseiny Sohi et al., 2014), - Estimation of permeability of grout curtain performed at the dam abutments and foundation by means of Lugeon tests- and cement take results (Hosseiny Sohi et al.,2017a), and - Monitoring the ground water level changes by means of standpipe piezometers embedded on the grouting galleries and drilled observation wells located in different places around the dam area (Hosseiny Sohi et al., 2017b). The findings and conclusions from these studies provide us insight and corroborate the need for developing proper numerical models to simulate ground water flow in order to: - Study the situation and patterns of the groundwater levels in the area, - Inspect the existing faults as possible water paths; and, - Study different scenarios regarding the water leakage. In the present study, using the topography- and geological maps of the dam site, reservoir water levels, Lugeon test data, piezometer-, observation well- and drainage borehole- data, a transient numerical model will be designed, constructed, calibrated and validated. The developed model simulates the dual porous media of Karun 4 dam using MODFLOW-2005 including Conduit 53 Flow Process (CFP) applied to the matrix medium and conduits, respectively. The process of calibration is performed in parallel for pipe calibration parameters and matrix medium. After successful calibration and validation of the model, the permeability of the grout curtain is evaluated and the adjusted pipes results are presented. This model could then be employed to study different scenarios regarding the water leakage through the existing faults and hypothetical paths, which is the subject of a subsequent study of the authors. 4.2 Study area 4.2.1 Site specifications and geological setting The Karun 4 dam is located at geographical latitude 50° 28´ 19´´and longitude 31° 35´ 50´´ close to the source of the Karun River in southwest Iran, where other important national dams constructed on the Karun river are under operation. The dam reservoir has an area of 29 km2 area and a volume that reaches 2.2 billion cubic meters, since the maximum normal water level is 1028 m.a.s.l. The dam was constructed to regulate the Karun river’s discharge and to control the frequent destructive floods of the Karun River, but also to generate about 1GW hydroelectric power (Figure 4.1). The Karun 4 dam is located in the tectonically folded belt of the Zagros Mountains within a canyon that is situated in the southwestern flank of the Kuh-Sefid anticline which has a NW- Figure 4.1. Plan view of the Karun 4 dam and grout curtain in Iran 54 SE trend. Most of the reservoir lies on Cretaceous to Miocene limestone and marly limestone (Haghnejad et al., 2014). The geological formations around the dam site consist of carbonate layers of the Asmari formation (As) which can be divided into lower (As1), middle (As2) and upper (As3) units and marlstone and marly limestone of impervious Pabdeh formation (Pd) upstream of the dam site (see Figure 4.1). The geology of the dam site is affected by many faults. The lower unit of the Asmari formation (AS1) consists of limestone, marly limestone and marlstone with some porous limestone (Mahab Ghods Report, 2002), and can be subdivided from upstream to downstream into 10 more detailed subunits AS1-a to AS1-j. The lithological description of these subunits listed in Table 4.1. Table 4.1. Lithological specifications of Asmari 1 subunits Name Description AS1-a Limestone (about 90%) gray to brownish gray, thick to very thick, strength to very strength, with no considerable karstic phenomena. There is some surficial erosion which is related to marly limestone (about10%) AS1-b Limy marlstone mostly covered by colluviums and slopewash. This rock is gray to dark green in natural color, weak to relatively strength and with no karstic features AS1-c Limestone, thick to very thick bed, gray to light brown with no particular karstic phenomena. AS1-d More or less like AS1-b AS1-e Limestone with some (in unweathered locations) gray to blue outcrops (marly limestone) very thick layers and without considerable karstic phenomena AS1-f This unit commences from As1-e and continues to the first porous key bed. Composed mostly with limestone having slight to medium karstic features (vug). There are some (10%) porous limestones. AS1-g Porous limestone which are solution vugs gray to light brownish gray in natural colour. At the right abutment the fresh outcrop is gray. Regarding with porous characteristics, it is considered as key bed. AS1-h Limestone with some grey to blue outcrops (probably marly limestone very thickbed. There are more or less karstic vugs at the top of this sub-unit AS1-i This sub-unit is more or less similar to AS1-g so it is considered as key bed. AS1-j Lime stone with marly limestone. the Marl limestone is frequent in the lower part, thick to very thick, no important karstic features is seen in this sub-unit 4.2.2 Seepage controlling system To seal the dam foundation and abutments, a very large grout curtain was constructed. This construction was done through the excavation of 5 series of galleries in both left and right abutments at different elevations from the dam crest elevation to the bottom and, finally, the 55 Figure 4.2. Grout curtain expanded cross section foundation gallery, called the 806-gallery (see Figure 4.2). The boreholes were drilled, examined and injected with grout from inside the galleries to create an integrated barrier wall or, called grout curtain, against water leakage (Houlsby, 1990). Some of the primary drilled holes were selected as exploratory boreholes and the rock permeability tests and geological logs were carried out within the boreholes. They were then considered in the grouting phase as normal drilled boreholes and the amounts of their cement injected recorded. The data obtained from the exploratory boreholes were utilized to design and construct the grout curtain in two lines in order to sew/stitch the dam structure to the upstream Pabdeh impermeable formation. The seepage controlling system consists of the grout- and the drainage curtain with the latter made up of boreholes drilled at 40 meters’ horizontal distance downstream in two galleries in the left and right abutments and allows not only to detect and to explore the still leaking water through the grout curtain, but also aids to reduce the uplift pressure and so prevents the breaking of mass rock downstream. To monitor the dam’s hydraulic behavior, a network of instruments is installed in the dam body, grouting galleries and the different stations of the dam site. The most important instruments of the seepage controlling system are standpipe piezometers and observation wells. The former, also known as Casagrande piezometers, are installed in the grouting galleries upstream and downstream of the grout curtain, whereas the observation wells are drilled in the different points of the dam body, abutments, the drainage curtain and downstream of the dam site. 56 4.3 Numerical model 4.3.1 Conduit Flow Process (CFP) The applied mathematical model used in the present study is Conduit Flow Process (CFP) imbedded in the modular finite-difference ground-water flow model MODFLOW-2005. The model has been developed by Shoemaker et al (2008) from U.S. Geological Survey. The three modes of its model-use described by Shoemaker et al (2008) are: “CFP has the ability to simulate turbulent or laminar groundwater flow conditions by: 1- coupling the traditional groundwater flow equation with formulations for a 1- dimensional discrete network of cylindrical pipes (Mode 1, CFPM1), 2- inserting a high-conductivity flow layer that can switch between laminar and turbulent flow (Mode 2, CFPM2), and 3- simultaneously coupling a discrete pipe network while inserting a high-conductivity flow layer that can switch between laminar and turbulent flow (Mode 3, CFPM3).” As the developers pointed out, CFPM1 accounts for dissolution or biological burrowing features in carbonate aquifers and voids in fractured rock. Thus, Mode 1 is employed to simulate the karstificated limestone formation of Asmari 1 in Karun 4 dam foundation and abutments. 4.3.2 Governing equation for porous medium The governing equation applied to the porous medium is the general groundwater equation derived from Darcy law (1856) and the mass-conservation law as following: The Darcy’s law for a one-dimensional flow in a prism of porous material (Figure 4.3) is defined as (McDonald and Harbaugh, 1988): 𝑄 = 𝐾𝐴(ℎ1−ℎ2) 𝐿 (4.1) where K [LT-1] is the hydraulic conductivity of the matrix in the direction of flow A [L2] is the cross-sectional area perpendicular to the flow ℎ1 − ℎ2 [L] is the head difference across the prism parallel to flow L [L] is the length of the prism parallel to the flow path. 57 Figure 4.3 Darcy’s law for a prism of porous material (from McDonald and Harbaugh, 1988) From the continuity equation could be deduced that the sum of all flows into and out of the cell must be equal to the rate of change in storage within the cell. Under the assumption that the density of ground water is constant, the continuity equation expressing the balance of flow for a cell is (McDonald and Harbaugh, 1988): ∑ 𝑄𝑖 = 𝑆𝑠 𝜕ℎ 𝜕𝑡 ∆𝑉 (4.2) where Qi [L3T-1] is the flow rate across face i of the cell (= positive for inflow and negative for outflow) Ss [L-1] is specific storage, ΔV [L3] is the volume of the cell, and ℎ [L] is the potentiometric head; 𝑡 [T] is time. The term on the right-hand side is equivalent to the net change of volume of water taken into/from the storage over a time interval 𝜕t given a change in head of 𝜕h. Combining equations (4.1) and (4.2) together, the general governing differential equation can be derived (Anderson et al., 2015), which represent three-dimensional transient groundwater flow for heterogeneous and anisotropic conditions (McDonald and Harbaugh, 1988): 𝜕 𝜕𝑥 (𝐾𝑥𝑥 𝜕ℎ 𝜕𝑥 ) + 𝜕 𝜕𝑦 (𝐾𝑦𝑦 𝜕ℎ 𝜕𝑦 ) + 𝜕 𝜕𝑧 (𝐾𝑧𝑧 𝜕ℎ 𝜕𝑧 ) ± 𝑊 = 𝑆𝑠 ( 𝜕ℎ 𝜕𝑡 ) (4.3) where: 𝐾𝑥𝑥, 𝐾𝑦𝑦, and 𝐾𝑧𝑧 [LT-1] are the values of hydraulic conductivity along the x, y, and z axes respectively; 58 ℎ [L] is the potentiometric head; 𝑊 [T-1] is a volumetric source flux per unit volume; 𝑆𝑠 [L -1] is specific storage. This equation is applied for the numerical modelling by MODFLOW packages and assumes that ground water is incompressible with constant density and viscosity. 4.3.3 Governing equations for pipe flow The Darcy-Weisbach equation represents the water head loss along a pipe, which is applicable to both turbulent and laminar flow: ∆ℎ = ℎ𝐿 = 𝑓 ∆𝑙 𝑑 𝑉2 2𝑔 (4.4) where: ∆ℎ or ℎ𝐿 [L] is the head loss measured along the pipe of length ∆l [L], 𝑓 [unitless] is the friction factor, 𝑑 [L] is pipe diameter, 𝑉 [LT-1] is the mean velocity, and 𝑔 [LT-2] is the gravitational acceleration constant. Laminar frictional pipe flow can be approximated with the Hagen-Poiseuille equation (Vennard and Street 1975). Because the mean pipe flow velocity, V, is equal to the volumetric flow, Q [L3T-1], divided by the cross-sectional area perpendicular to flow, A [L2], the Hagen-Poiseuille equation can be reformulated to solve for volumetric flow rates, Q, rather than head losses Δh; was written in terms of the head gradient along the pipe, the viscosity of the fluid, and the pipe diameter squared by Craft and others (1962) and Halford (2000) as follows: 𝑄 = −A 𝜌𝑔𝑑2∆ℎ 32𝜇∆𝑙𝜏 (4.5) where 𝜇 [ML-3] is the density of the water; 𝜌 [ML-1T-1] is the dynamic viscosity of water; where τ is the tortuosity [unitless] of the pipe, which is defined as the distance actually traveled by a particle flowing through the center of the pipe divided by the straight line distance between the nodes of the pipe. 59 The applied relation between f and Reynolds number (Re) in CFPM1 for turbulent flow derived from Colebrook-White formula (Colebrook and White, 1937; Dreybrodt, 1988): 1 √𝑓 = −2 log( 𝐾𝑐 3.71𝑑 + 2.51 𝑅𝑒√𝑓 ) (4.6) where: 𝑓 [unitless] is the friction factor; 𝐾𝑐 [L] is the mean roughness height of the conduit wall; and 𝑑 [L] is pipe diameter. 𝑅𝑒 is defined as: 𝑅𝑒 = 𝑉𝑑𝑝𝑜𝑟𝑒𝜌 𝜇 = 𝑉𝑑𝑝𝑜𝑟𝑒 𝜈 (4.7) where: 𝑑𝑝𝑜𝑟𝑒 [L] is the mean pore size diameter, and, 𝜈 [L2T-1] is the kinematic viscosity of water. The exchange discharge between matrix medium represented by MODFLOW and conduits is computed through the following linear equation in CFP (Shoemaker et al.,2008): 𝑄𝑒𝑥 = 𝛼𝑗,𝑖,𝑘(ℎ𝑖𝑛 − ℎ𝑗,𝑖,𝑘) (4.8) Where: 𝑄𝑒𝑥 [L3T-1] is the volumetric exchange flow rate, 𝛼𝑗,𝑖,𝑘 [L 2T-1] is the pipe conductance at MODFLOW cell j,i,k ℎ𝑖𝑛 [L] is the head at the pipe node n located at the center of the MODFLOW cell, ℎ𝑗,𝑖,𝑘 [L] is the head in the encompassing MODFLOW cell j,i,k, and 𝛼𝑗,𝑖,𝑘 is computed by CFP as: 𝛼𝑗,𝑖,𝑘 = ∑ 𝐾𝑗,𝑖,𝑘𝜋𝑑𝑖𝑝 1 2 (Δ𝑙𝑖𝑝𝜏𝑖𝑝) 𝑛𝑝 𝑖𝑝=1 𝑟𝑖𝑝 (4.9) where 𝐾𝑗,𝑖,𝑘 [LT -1] is the conduit wall permeability (conductivity) in MODFLOW cell j,i,k, 𝑑𝑖𝑝 [L] is the diameter of pipe ip connected to node in, 𝛥𝑙𝑖𝑝 [L] is the straight-line length between nodes of pipe ip connected to node in, 60 𝜏𝑖𝑝 [unitless] is the tortuosity of pipe ip connected to node in, and 𝑟𝑖𝑝 [L] is the radius of the pipe [L]. 4.3.4 Graphical User Interface ModelMuse The graphical user interface ModelMuse for MODFLOW-2005as developed by Richard B. Winston (2009) of USGS, version 3, which supports CFP has been applied to set up the Karun 4 dam numerical model. 4.3.5 Set-up and designing of the numerical model 4.3.5.1 Conceptual model Figure 4.4 presents the conceptual model of the dam area and the most possible groundwater flow paths. The general paths of the groundwater flow from the reservoir downstream the river valley are illustrated. Also shown are the boundary conditions to be discussed later. The conceptual model has been developed based on the collected site data, and site specifications. 4.3.5.2 Model grid To build the model grid, the drawing (DWG format) received from the Iran Water and Power Resources Development Company (IWPC) are imported into the ARC-GIS environment to provide the corresponding shape- or ASCII files. The ModelMuse is used to assign the layout and attributes of the civil structures of the dam area, such as, dam body, grouting curtain, drainage curtain, observation wells and piezometers. Starting with a rather coarse spatial discretization of the conceptual model with a 25x25 m2 cell size, this grid has then been refined to 5x5 m2 in subsequent steps. The entire study area includes 170 rows and 150 columns in north-south and west-east directions, respectively, i.e. a total of 25500 cells in the horizontal plan view. However, as listed in Table 4.2, the number of active cells in each layer is much less. The north side of the grid is generally limited by the impervious Pabdeh formation. Only a small part of the grid is dedicated to the reservoir’s contact with the topography of the valley. Both the west and east sides of the grid are restricted to the elevation of 1035m since the 61 Figure 4.4. Conceptual model of Karun 4 dam groundwater flow and boundary conditions layout maximum normal water level of the dam is 1028m. The dam crest elevation itself has been set at an elevation of 1032m. The south side where the line connects the furthest observation wells at the left and right banks is modeled with an offset of 50 meters. These limitations are used to avoid having unnecessary cells that would cause extreme computational time of the model, particularly, in transient state. As mentioned in Section 4.2.1, the geological characteristics of the Karun dam site are too diverse, and change drastically from one subunit of Asmari to the neighboring site. Hence, only one layer could not represent the real circumstances in the vertical direction. To account for this site variation characteristics, four layers are defined, as indicated in Table 4.2, based on the elevation of the excavated grouting galleries, where the standpipe piezometers are installed for measurements purposes. The topography is assigned to the top layer. Table 4.2. Layer specifications assigned to the aquifer model Name Elevation (m.a.s.l) Type of Aquifer Number of active cells Bottom Top Layer 1 (top layer) 974 1032/topography unconfined 6,246 Layer 2 890 974 confined 7,747 Layer 3 848 890 confined 9,248 Layer 4 (bottom layer) 700/Pabdeh 848 Confined 10,217 62 Figure 4.5. Plan view of the 4 layers of the Karun 4 MODFLOW discretized model with the corresponding bottom and top elevations. According to the monitoring reports from the dam site (Mahab Ghods Consulting Engineering, 2015), the piezometers in the elevations lower than 974 meter show groundwater heads higher than their installed elevation (artesian pressure), therefore, the three bottom layers of model are defined as confined layers, whereas the top layer as unconfined. The plan view of the model grid for the four layers is illustrated in the previously shown Figure 4.5. 63 Figure 4.6. Reservoir water level hydrograph from time of impounding till end of modelling time 4.3.5.3 Boundary and initial conditions Boundary conditions of model are defined as in the following and presented in Figure 4.4. Reservoir: CHD (Time-Variant Specified-Head Package), The reservoir water level is measured daily at the dam site and is utilized as input to the model. This measurement is the only independent variable that changes during time steps and affects the entire model. Figure 4.6 shows the time-series of the reservoir water level from first days of impounding in March 2010 till March 2017 (end of modelling duration). River: CHD (Time-Variant Specified-Head Package), The Karun river is defined as CHD boundary with a constant head of 845 m.a.s.l. Pabdeh impermeable formation borders in upstream: no flow (inactive cells), Downstream and lateral borders: CHD (derived from the extrapolated head measurements). After calibrating the steady state phase, the output of simulated values of groundwater head have been assigned as initial conditions of groundwater head in the transient state. 4.3.5.4 Hydraulic conductivity The geological formations of the model area include Pabdeh- and Asmari1 formation; and, Asmari 2 formation in some areas at downstream of the reservoir. The lower Asmari formation (AS1) subunits are defined as the local zones in each layer to calibrate the model by means of the hydraulic conductivity (from now onwards: K). The zones have been extracted 64 Figure 4.7. Geological map for layer 2 including the Asmari 1 subunits from the geological maps of IWPC. The initial values of K are obtained from the hydraulic conductivity (see Table 4.3) proposed by the consulting engineers of the Karun 4 project (Sadd Tunnel Pars Company Report, 2015). Figure 4.7 displays, the geological map of layer 2 that is employed to assign the K zones. Table 4.3. Proposed hydraulic conductivity values by Sadd Tunnel Pars company (2015) Description K (m/s) K (m/day) Grout curtain 9.1 x 10-7 0.079 Drainage curtain 1.3 x 10-4 11.232 Pabdeh 6.5 x 10-7 0.056 AS1-a, AS1-c, AS1-e, AS1-f, AS1-j 2.6 x 10-5 2.246 AS1-b1 1.3 x 10-4 11.232 AS1-b 2.6 x 10-6 0.225 AS1-d 1.3 x 10-6 0.112 AS1-g, AS1-i 1.3 x 10-4 11.232 AS1-h 5.2 x 10-6 0.449 4.3.5.5 Drainage conductance The second calibration parameter is the hydraulic conductance, C of the drainage curtain using the drainage package (DRN) in MODFLOW-2005. A drain (Figure 4.8) removes water from the aquifer at a rate proportional to the difference between the head in the aquifer and 65 Figure 4.8 Buried drain pipe (from McDonald and Harbaugh, 1988) the drain elevation, as long as the head in the aquifer is above that elevation. Based on this concept the drain conductance CD is defined over the following equation (Harbaugh, 2005): Q out= CD (hi,j,k – HD), when hi,j,k > HD (4.10a) and Q out = 0 when hi,j,k ≤ HD (4.10b) where Q out [L3T-1] is the flow from the aquifer into the drain, CD [L2T-1] is the drain conductance HD [L] is the drain elevation at the borehole bottom, and hi,j,k [ L] is the head in the cell containing the drain. The conductance CD can be, based on Darcy’s law, composed further as CD = 𝐾𝐴 𝐿 (4.11) where K [LT-1] , A [L2] and L [L] are the hydraulic conductivity, mantle area and the length of the pipe, respectively. Because of this equation, CD can be calibrated over the pipe’s hydraulic conductivity. 66 Figure 4.9. The calibrated zoning sections of the drainage curtain To calibrate the model more precisely, the drainage curtain is divided into 10 different sections/sectors, as shown in Figure 4.9. The drainage elevations at the beginning sector in both banks close to the river are assigned 762 m.a.s.l and the remaining sectors 812 m.a.s.l. 4.3.5.6 Head observation package (HOB) The head observation package (HOB) is used to compare simulated and observed heads in the observation wells and stand-pipe piezometers, with their locations provided in Table 4.4. Table 4.4. Observation wells’ and piezometers’ specifications 67 Figure 4.10 Layout of piezometers and observation wells Figure 4.10 shows the layout of the observation wells and piezometers in the plan view. 4.3.5.7 Storage Specific storage and specific yield are to be calibrated for the different layers in the model in the transient phase. Specific storage and specific yield make together the storativity. In a confined aquifer, the storativity is defined as (Anderson et al., 2015): S = Ss b (4.12a) where Ss [L-1] is the specific storage, and b [L] is the aquifer/layer thickness. In an unconfined aquifer, the storativity I is defined as: S = Sy+ Ss b ≅ Sy (4.12b) where Sy [unitless] is the specific yield. Since in an unconfined aquifer Sy is much larger than Ss b, its storativity is effectively equal to Sy . The zonation of the storativity is the same as that of K and corresponds to the formations of the aforementioned geological maps of the dam site. 68 4.4 Model calibration 4.4.1 Calibration parameters of porous medium in steady state For steady-state conditions, time plays no role on the outcome of the model. Hence the pipes are not modeled in this phase, which makes the numerical groundwater flow model a common MODFLOW model. The steady-state calibration is achieved by using the collected data, the assumptions described in the previous sections, and by adjusting the two calibration parameters K and C through a trial and error process, until the differences (residuals) between measured and simulated groundwater heads are somewhat small enough. The criterion for the goodness of calibration is the root mean squared residual (rmsr) calculated as: rmsr= √ ∑ (𝑅𝑖)2 𝑛 1 𝑛 (4.13) where Ri is the residual between observed- and simulated head datum i in the piezometer/observation well and n is the total number of observed data, i.e. the product of the number of piezometers/observation wells (=23) and the number of stress (observation) periods which is =1 for the steady-state calibration – for which the heads on 20.06.2012 are taken as the calibration targets - and =59 for the transient calibration. The zonal calibration of K, particularly, close to the grout curtain wings, downstream of the reservoir in layer 3 and 4, is performed not only by changing the absolute values, but also by further sub-division, i.e. refinement of these zones. 4.4.2 Calibration parameters of porous medium in transient state Assigning the steady- state-calibrated K and C, the storativities Ss and Sy are calibrated in the local zones during the transient calibration process by a trial and error. Unlike for steady state, the head/flow conditions and specifications change during the transient state. As mentioned the head data observed at the dam site covers an interval of the 57 months, which are divided into two parts for calibration and validation, respectively. The start time of the calibration period is June 20, 2012 and includes 38 stress periods, up to June 18, 2015. The remaining 21 stress periods cover then the validation period, up to March 13, 2017. Table 4.5 lists the stress-periods/dates/lengths used. 69 Table 4.5. Stress periods with dates and lengths of transient model 4.4.3 CFP calibration (transient state) After calibration of the entire, continuum porous medium- MODFLOW model, the potential conduits/faults are put into the model at the locations of the recognized faults and the CFP- calibration process is started. Figure 4.11 shows the layout of faults in the dam plan view. Since the conduits of the karstificated Asmari 1 formation can be assumed as pipes, the conduit flow pipes mode (CFPM1) is chosen to simulate the conduits of the dual porosity medium of the Karun 4 dam. The layout of the pipes is extracted from the Karun 4 dam geological maps (IWPC) (see Figure 4.11). These drawings are prepared at the corresponding elevations of the galleries of the grout curtain, which is in harmony with the defined layers of the present model. Stress Period Date Day number Step length (day) Stress Period Date Day number Step length (day) 1 20.06.2012 0 Steady state 39 17.07.2015 1122 29 2 26.07.2012 36 36 40 19.08.2015 1155 33 3 20.08.2012 61 25 41 18.09.2015 1185 30 4 17.09.2012 89 28 42 23.10.2015 1220 35 5 26.10.2012 128 39 43 20.11.2015 1248 28 6 26.11.2012 159 31 44 18.12.2015 1276 28 7 24.12.2012 187 28 45 15.01.2016 1304 28 8 21.01.2013 215 28 46 12.02.2016 1332 28 9 04.03.2013 257 42 47 12.03.2016 1361 29 10 01.04.2013 285 28 48 17.04.2016 1397 36 11 29.04.2013 313 28 49 14.05.2016 1424 27 12 28.05.2013 342 29 50 18.06.2016 1459 35 13 01.07.2013 376 34 51 16.07.2016 1487 28 14 06.08.2013 412 36 52 12.08.2016 1514 27 15 07.09.2013 444 32 53 09.09.2016 1542 28 16 07.10.2013 474 30 54 14.10.2016 1577 35 17 12.11.2013 510 36 55 19.11.2016 1613 36 18 10.12.2013 538 28 56 16.12.2016 1640 27 19 05.01.2014 564 26 57 16.01.2017 1671 31 20 29.01.2014 588 24 58 13.02.2017 1699 28 21 19.02.2014 609 21 59 13.03.2017 1727 28 22 15.03.2014 633 24 23 12.04.2014 661 28 24 18.05.2014 697 36 25 05.06.2014 715 18 26 11.07.2014 751 36 27 10.08.2014 781 30 28 04.09.2014 806 25 29 09.10.2014 841 35 30 06.11.2014 869 28 31 11.12.2014 904 35 32 04.01.2015 928 24 33 22.01.2015 946 18 34 19.02.2015 974 28 35 22.03.2015 1005 31 36 23.04.2015 1037 32 37 21.05.2015 1065 28 38 18.06.2015 1093 28 ValidationCalibration 70 Figure 4.11 Karun4 dam area plan view with faults sketched in red According to Eq. (4.9), the corresponding CFP parameters of each pipe to assign are the followings: pipe diameter, tortuosity, roughness height, lower- and higher critical Reynolds number, conduit wall conductivity term and elevation of the pipe. However, the calibration parameters of the pipes are only its diameter and conduit wall conductivity. The former enters in the various equations described in Section 4.3.3, i.e. the Reynolds number, determining the flow regime in the pipe. Hence, if the flow is laminar, the Hagen-Poiseuille equation (4.5) is solved, whereas for turbulent flow, the Colebrook-White formula (4.6) will be applied to solve the Darcy-Weisbach equation (4.4). The conduit wall conductivity determines the exchange discharge between the conduit and the matrix medium and enters Eqs. (4.8) and (4.9). Other parameters of Eq. (4.9) have no significant impact on the calibration process. Thus, the following values are assigned to all pipes: Tortuosity: 1.0 Roughness height: 0.01 Lower critical Reynolds number: 2 Higher critical Reynolds number: 10 (Shoemaker et al., 2008) Elevation of pipe: 895, 850 and 800 m.a.s.l to layers 2, 3 and 4, respectively. 71 The calibration is done for each pipe to find out the appropriate diameter and conductance considering the above-mentioned assumptions. Firstly, the diameter of a pipe, starting from D=0.001 meter, is determined , by running the model and computing its rmsr. If the latter is reducing, D moves into the right direction. This process of increasing D is continued, until the rmsr of the model does not decrease anymore. Subsequently, the conductivity of the pipe is calibrated in a similar way, starting from K= 0.00001 up to 0.12 m/day. Here again, the proper conductance should reduce the rmsr. When this step is finished, the fault is adjusted to the model as a pipe; and the same process will be proceeded for the other faults. In both the pipe diameter D- and K- calibration, the only criteria to decide to keep the corresponding value, or even the fault as a pipe, was its effect on the rmsr of the model, i.e. reducing it. Therefore, as not all candidate faults retrieved from the geological maps could be calibrated as discussed above, they were discarded from the CFP- model. Consequently, the finally defined 15 pipes of the model illustrated in Figure 4.12 describe only those faults whose characteristics can be most likely expressed as a discrete conduit within the matrix medium in a karstified formation to be simulated with the CFPM1- model mode. At the left bank there are more matching pipes, because the drainage there is almost three times higher and the variation of the observation wells’ residuals is higher. This confirms the results of the previous study of the authors (Hosseiny Sohi et al, 2017a). In the upper layers, especially layer 2, toward downstream of the reservoir, some of the pipes intend to be dry. It should be avoided to have the dry pipes, which causes errors in the model, by limiting the length of the critical pipe to the wet area. The wet area can be obtained from the water head contours which is available in the results section by comparing their numbers with those of the top and bottom elevations of the corresponding layer. 72 Figure 4.12. Layout of the pipes at a) layer 2, b) layer 3, and c) layer 4, derived from the locations of the detected faults based on the geological maps, d) faults plan view 4.5 Results 4.5.1 Steady state calibration 4.5.1.1 Hydraulic conductivity (K) Using the K-zoning described in Section 4.3.5.4, the finally calibrated K-values in steady-state mode are shown for the four layers in Figure 4.13. 73 4.5.1.2 Drain conductance (C) The drainage conductance C (see Eq. 4.10) for each sector is calibrated separately using as calibration target the total observed values of the drained water discharge at the left- and right drainage galleries and which are equal to 0.4 and 0.13 m3/s, respectively (Mahab Ghods monitoring report, 2015). After several try-and-error calibration runs of the MODFLOW- model, the calibrated values of C as listed in Table 4.6 are obtained. Table 4.6 Characteristics of the different sections of the drainage curtain. Drain ID Length (m) Conductance (C) (m2/day) Drain_1 63.7 34.6 Drain_2 55.0 14.5 Drain_3 72.6 2.59 Drain_4 78.3 0.09 Drain_5 74.2 1.73 Drain_6 72.3 0.86 Drain_7 80.6 1.56 Drain_8 39.8 1.04 Drain_9 58.4 2.76 Drain_10 133.8 3.89 Different geological characteristics of Asmari 1 subunits can explain the big difference among the diverse sections of drainage curtain. The sections are located at the intersection of AS1-g and AS1-i show more water drainage, which is a consequence of the high permeability of these areas. Furthermore, the conductance on the right side is generally lower than on the left side, which is again a confirmation of the geological specification of the region. 4.5.1.3 Groundwater head The rmsr of the calibrated steady-state model of groundwater head data reaches 1.44, with the highest residuals remaining in the range of ±2.21m. Figure 4.14 (left panel) shows the simulated- over the observed heads in the calibration period and one can notice that these lie very well on the slope =1 regression line, as quantified by the very high coefficient of determination, or squared correlation coefficient of 0.998 between the two data sets. The residuals against the observed values are illustrated in Figure 4.14 (right panel). The isolines of the simulated groundwater heads in layer 3 are demonstrated in Figure 4.15. One can see that their pattern has a close similarity with those of the earlier conceptual model, which attests to some degree the validity of the numerical model. 74 Figure 4.13 Calibrated K-values of in a) layer 1, b) layer 2, c) layer 3, and, d) layer 4. 75 Figure 4.14. Steady-state calibrated- versus observed heads (left panel) and residuals (right panel) Figure 4.15. Steady-state calibrated ground water heads in layer 3 4.5.2 Transient state calibration 4.5.2.1 Storage As the top layer is an unconfined layer, the highest value of specific yield (Sy) equal to 0.1 is assigned to this layer which improves the rmsr of the model. For the three confined bottom layer Sy =1E-7 is specified, as they have essentially no specific yield. The specific storage Ss of all zones of the top layer is specified as 0.02 m-1, because changing Ss in different zones of this layer does not have a significant impact on the model outcome. The calibrated values Ss for layers 2, 3 and 4 are shown in Figure 4.16. 76 Legend Figure 4.16 Calibrated Ss-values of in the (confined) a) layer 2, b) layer 3, c) layer 4. 4.5.2.2 Grout-curtain-efficiency evaluation The engineers of Sad Tunnel Pars company (2015) proposed one constant value of K= 9.1E-7 m/s, equal to 0,079 m/day, for the hydraulic conductivity (or tightness) along the entire grout curtain. From the calibration process and having access to observed measurements in the piezometers and wells, a variable distribution of K in each layer and dam wing is obtained here, instead. Therefore, both left and right wings in each layer are divided into two parts, with the 77 first one (AB and DE) connected to the dam body and the second one (BC and EF) running more or less perpendicular to it in upstream direction toward the Pabdeh formation (see Figure 4.18). Table 4.7 lists the finally calibrated K-values in these two parts and one can notice that almost all of them are higher than those assumed by the design engineers, particularly, in the left wing. Thus, the advantage of the use of the calibrated MODFLOW-CFP numerical model here is to be able to better evaluate the practical uncertainties in the construction of the grout curtain, as manifested by the Lugeon tests and cement Take, as described in detail in that earlier study of the authors (Hosseiny Sohi et al., 2017a). Table 4.7. Hydraulic conductivity (m/day) of the grout curtain derived from the calibrate model Layer Left wing connected to: Right wing connected to: Dam body(AB) Pabdeh formation(BC) Dam body(DE) Pabdeh formation(EF) 1 0.092 0.092 0.092 0.084 2 0.092 0.092 0.092 0.084 3 0.101 0.092 0.075 0.084 4 0.092 0.075 0.086 0.084 4.5.2.3 Groundwater heads Applying all above-mentioned steps, the transient calibration process has been performed. In this case the whole time series from June 2012 to March 2017 has been divided in one set for the calibration from June 2012 to June 2015 and the remainder for the validation (June 2015 to March 2017). The simulated and observed groundwater heads are compared by means of the head observation package (HOB) of MODFLOW-2005. Figure 4.17 shows the results for some of the selected piezometers for both the calibration- and validation phase. One can notice that although the simulated and observed head time series are running very much in parallel, depending on the locations of the piezometer or observation wells, there are systematic over- or underestimations of the observed head curves, i.e. the residuals are consistently negative or positive, respectively. To understand if these head-residual trends are somehow related to the location of the observation point within the Asmari 1 formation, the most extreme station residuals are displayed as red (negative residuals) and blue (positive ones) in Figure 4.18. Obviously, no spatial trend is visible, as the head residuals change abruptly across neighboring observations locations, most likely due to strong geological heterogeneities in the fractured karst dam area. 78 G ro u n d w a te r h ea d ( m .a .s .l ) Figure 4.17 Observed- and simulated groundwater head time series in the calibration phase form June 2012 to June 2015(left panels) and in the validation phase from June 2015 to March 2017 ( right panels) for different observation points: a) OW1, b) OW4, c) OW6 and d) P4 79 Figure 4.18 Residual of observed- and simulated groundwater head (m) of all piezometers and observation wells (negative ones in red), (positive ones in blue) Figure 4.19 shows the observed- versus simulated values (left panel) and the residuals against the observed measurements (right panel) for the calibration phase (June 2012 to June 2015) and Figure 4.20 for the validation phase. The R2 of the slope=1 regression line of the simulated- versus observed groundwater heads reaches 0.98 and 0.97 in calibration and validation phase, respectively. The results indicate a very strong correlation between simulated- and observed head values for both calibration- and validation phase, wherefore, as is to be expected, R2 is higher for the former than for the latter. When taking all 59 stress periods for the calibration and the validation together, the rmsr of whole final model is rmsr =2.44 m and all residuals < ± 6.80m, and these are the reference values for the evaluation of any further model scenario, as discussed in the next chapter. The simulated- versus observed groundwater head in OW1 (left) and OW4(right) as typical results of observation wells are displayed in Figure 4.21. The coefficient of determination (R2) of OW1 data is 0.98 and that of OW4 is 0.91. 80 Figure 4.19 Simulated- versus observed head values (left) and residuals of each observation point (right) – Transient state (calibration phase) Figure 4.20 Similar to Figure 4.19, but for the validation phase Figure 4.21 Similar to Figure 4.19 (left), but for well OW1 and OW4 for the calibration phase 81 Figure 4.22 Simulated groundwater head contours of the last stress period on March 13th 2017 in : a) layer 2, b) layer 3 and c) layer 4. Figure 4.22 shows the groundwater flow plan view of the layer 2, 3 and 4 on the 59th (last) stress period (March 13th 2017). The water in the unconfined layer 1 exists only in the reservoir,and in its dowanstream areas all cells are dry. As shown in this figure, the groundwater level in the reservoir is very high, 1028 (maxium reservoir water level) which decreases to 845 m.s.a.l (river water level) toward the downstream boundaries. The intense contours around the grout curtain are a consequence of its function. In layer 4, toward downstream, the impact of the drainage curtian is clearly traceable. 82 4.5.3 Impact of faults on the groundwater leakage through the matrix The faults and fractures are well-known possible water leakage paths through the limestone and carbonate formation with its karstification potential. The implemented grout curtain would not necessarily seal the water leakage entirely and the faults remain as water paths through the dam area. To inspect whether a fault passes water flow or not, a pipe is set at its location. The variable calibration parameters of the pipe are diameter and conduit wall permeability (conductivity), with their final values listed in Table 4.8. Each pipe/fault includes several corresponding nodes and tubes in the model that are produced form the intersected cells of the grid. For example, pipe no. 7 from layer 3, contains 39 nodes and 138 meters length. Table 4.8. Pipe specifications assigned to the CFP- model (see Figure 4.12 for locations) The various panels of Figure 4.23 illustrate the water heads in both nodes and the matrix medium along the pipe nodes for each of the 15 pipes (left y-axis) (see Table 4.8 and Figure 4.12 for locations), as well as the exchange flow Q between porous medium and pipe resulting from these head differences (right y-axis). Positive values of Q mean flow from the pipe to the matrix and vice versa (see equation 4.8). The fault F4-Up in the left abutment is represented in layer 4 by pipe 1 inside- and pipe 5 outside the reservoir; in layer 3 by pipe 8 inside- and pipe 9 outside the reservoir. Pipe 3 and pipe 4 are the two sides of fault F13 on the right abutment at the upstream and downstream of the reservoir, Elevation (m.a.s.l) Pipe K (m/day) Diameter (m) Length (m) Begin Node Final Node Fault Location 1 0.005 0.02 48 1 11 F4-Up Upstream 2 0.005 0.1 211 12 71 F7 Downstream 3 0.00005 0.01 116 72 91 F13 Upstream 4 0.09 0.005 106 92 113 F13 Downstream 5 0.1 0.02 310 114 184 F4-Up Downstream 6 0.1 0.1 93 185 207 F2 Upstream 7 0.12 0.079 138 208 246 F3 Downstream 8 0.01 0.007 47 247 257 F4-Up Upstream 9 0.05 0.01 299 258 325 F4-Up Downstream 10 0.09 0.05 283 326 389 F4-Down Downstream 11 0.0005 0.01 49 390 401 F2 Downstream 12 0.00001 0.05 43 402 413 F2 Upstream 13 0.00007 0.01 91 414 440 F3 Upstream 14 0.12 0.1 171 441 488 F15 Upstream 15 0.00001 0.05 101 489 517 F14 Upstream 848 (layer 4) 890 (layer 3) 932 (layer 2) 83 Figure 4.23 Water heads in the nodes of the various pipes and the porous matrix medium, and the exchange flow between them, on March 13th 2017, the last stress period of the model 84 Figure 4.23 (continued) 85 Figure 4.23 (continued) respectively, and pipes 12 and pipe 11 are the two sides of the fault F3 in layer 2 in the upstream and downstream of the reservoir, respectively. Therefore, the graphs of all these pair-pipes are depicted together in one panel. Only in five pipes are the node water heads proportional to those in the porous matrix: Pipe 5, pipe 8, pipe 9 and pipe 10 in the left abutment and pipe 4 in the right abutment. In all other pipes the water heads are almost constant, with changes < 0.5 m along their courses. 86 In layer 4, the pipes 4 and 5 are located downstream of the grout curtain on the right- and left abutments, respectively (Figure 4.12). Pipe 9 is located on top of pipe 5, but at elevation 890m (layer 4), whereas its diameter and K are half of those of pipe 5. Pipe 8, which is the continuation of pipe 9 toward upstream of the reservoir, has a smaller diameter and it’s K is only one fifth of that of pipe 9. This heads in this pipe are in significant harmony with those in the matrix. For pipe 10 in layer 3, close to pipe 9, but with a larger diameter D, and an even larger D than pipe 5 in the lower layer, and with a K factor two times bigger than pipe 9, the heads are again in tandem with those of the matrix. The pipes 1, 2, 3, 11, 12, 13 and 15 have a very low K- values, which means there is a very low exchange between matrix and pipe. The exchange flow value in these pipes is less than ± 2 m3/day. This is to identify through the corresponding graphs which the water head of the pipe is almost constant and do not follow the pattern of matrix at the same node. Pipe 4, despite its high K-value, has a low exchange flow which is explainable with its relatively low D-value. Pipe 6 (fault F2) and pipe 7 (fault F3) in layer 3 have not only large diameters but also high K- values. As the water head of these pipes are constant and is not similar to the matrix corresponding nodes, it could be concluded that these two faults F2 and F3 despite its high values of pipe diameter and K, have negligible impact on the water leakage, which also explains that this wing of the grout curtain (perpendicular to the Pabdeh formation toward upstream) seals water more than the wing is connected to the dam body. This behavior could be explained again by the position of these faults which are located close to impermeable Pabdeh formation. However, due to both high K- and D-values, the exchange flow in both pipes is high. Fault F4 Up, in contrast, represented by pipe 5 in layer 4, and pipe 8 and pipe 9 in layer 3 is the most critical fault in the left abutment and the entire model area. The highest fluctuation of the exchange flow between matrix and pipes happen there, even with low diameter (pipe 5) and low K (pipe 9). Pipe 10 displays the fault F4 Down, close to F4 Up at the left abutment (see Figure 4.12). The exchange flow in this fault changes from -15 to 5 m3/day. In the right abutment fault F15 represented by pipe 14, the heads of the matrix- and conduit nodes do not have the same head pattern again. Here the exchange flow between them reaches again high values, due to both high K- and D- values which are equally important to identify the amount of exchange flow between porous medium and conduits. 87 As there are many sophisticated factors which affect the water flow exchange between conduits and matrix medium, it is not simple to extract an explicit relation between them, although the results obtained and discussed above appear to indicate a high amount of matching with the statements in the geological reports of the designer engineers and the observed data from dam site. However, it is recommended to investigate and inspect in more detail the area close to fault F4 in layers 3 and 4 by the installation of new monitoring wells. 4.6 Conclusions The present model employs a novel approach in numerical modelling of groundwater flow in Karst aquifer, using MODFLOW-2005, the well-known code of USGS; in combination with a turbulent flow code named Conduit Flow Process CFP), in the site of a concrete double-arc dam built in a very narrow valley. Where the governing equation of MODFLOW is Darcy’s law, CFP uses Darcy-Weisbach general equation to solve the groundwater flow. The application of these codes to the Karun 4 dam project had various challenges which have been successfully met and solved. The simulated values of the groundwater flow of Karun 4 ground water have a significant correlation with observation measurements form the dam site in both steady- and transient state. In the steady state, a laminar groundwater model is set and calibrated in MODFLOW changing the hydraulic conductivity of the different zones of Asmari1 formation in four layers of the model; and, the conductance of the drainage curtain on the downstream of the reservoir. Then in the transient state, firstly, the matrix medium is calibrated in MODFLOW changing the storage parameters of the dam site. Afterwards, pipes are set to CFP and calibrated. The calibration parameters of the pipes are their hydraulic conductivity and diameter. From the recognized layout of the faults from the geological maps, the selected pipes are set to the model. The criterion of the selection is to reduce rmsr of the whole model. During the calibration process, the hydraulic conductivity values of the grout curtain are calibrated too, which indicated that the performed curtain is not a homogeneous obstacle against the water leakage. The faults located on the left abutment indicate a more significant impact on the water leakage, since the drainage amount are higher there. The most critical fault, namely fault F4-Up, is identified in this abutment and in lower elevations of the Asmari1 formation, which conducts 88 water from upstream of the reservoir toward downstream. It is recommended to monitor this fault more precisely by installing new observation wells around it. The different scenarios of water leakage could be studied using present calibrated model, which is the issue of an upcoming study. References Anderson MP, Woessner WW & Hunt RJ (2015). Applied groundwater modeling: Simulation of flow and advective transport, 2nd Edition. Assari A & Mohammadi Z (2017). Assessing flow paths in a karst aquifer based on multiple dye tracing tests using stochastic simulation and the MODFLOW-CFP code, Hydrogeology Journal, Vol 25, No 6, pp 1679-1702. Ashok KC & Sophocleous MA (2009). Recent MODFLOW developments for groundwater modeling: Kansas Geological Survey Open-file Report 2009-4, 139 p. Atkinson TC (1977). Diffuse flow and conduit flow in a limestone terrain in the Mendip Hills, Somerset, J. Hydro., 19, 323– 349. Barenblatt GK, Zheltov IP & Kochina N (1960). Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks, Prikl. Mat. Mekh., 24, 5, 852– 864. Clemens T, Huckinghaus D, Sauter M, Liedl R & Teusch G (1996). A combined continuum and discrete network reactive transport model for simulation of karst development. Proceedings of the ModelCARE 96 Conference held at Golden, Colorado, September 1996. IAHS Publ. no. 237, 306- 318. Colebrook CF & White CM (1937). The reduction of carrying capacity of pipes with age: London, Journal Institute of Civil Engineers, v. 7, p. 99. Cornaton F and Perrochet P (2002). Analytical 1-D dual-porosity equivalent solutions to 3-D discrete single continuum models. Application to karstic spring hydrograph modeling, J. Hydrology, 262, 165-176. Craft BC, Holden WR & Graves EDJr (1962). Well design: Drilling and production: Englewood Cliffs, N.J., Prentice-Hall, 571 p. Darcy H (1856). Les Fontances publiques de la ville de Dijon: Paris, Victor Dalmont. Doğan U & Çiçek I (2002). Occurrence of cover-collapse sinkholes in the May Dam reservoir area, Konya, Turkey, Cave and Karst Science, Vol. 29, No.3. Dreybrodt W (1988). Processes in karst systems—Physics, chemistry, and geology: Berlin, Springer- Verlag, 288 p. Ford D & Williams P (2007). Karst Hydrogeology and Geomorphology. Wiley Ltd, England. Gallegos JJ, Hu BX, & Davis H (2013). Simulating flow in karst aquifers at laboratory and sub-regional scales using MODFLOW-CFP: Hydrogeology Journal, v. 21, p. 1749-1760. 89 Ghasemizadeh R, Hellweger F, Butscher CH, Padilla I, Vesper D, Field M & Alshawabkeh A (2012). Review: Groundwater flow and transport modeling of karst aquifers, with particular reference to the North Coast Limestone aquifer system of Puerto Rico. Gunn J (1985). a conceptual model for conduit flow dominated karst aquifers, Karst Water Resources (Proceedings of the Ankara - Antalya Symposium, July 1985). IAHSPubl. no. 161. Haghnejad A, Ahangari K & Noorzad A (2014). Investigation on various relations between uniaxial compressive strength, elasticity and deformation modulus of Asmari formation in Iran. Arab. Journ. Science and Engineer., 39, 4, 2677-2682. Halford KJ (2000). Simulation and interpretation of borehole flowmeter results under laminar and turbulent flow conditions: Seventh International Symposium on Logging for Minerals and Geotechnical Applications, Golden, Colo., October 24–26, 2000, p. 157–168. Halihan T & Wicks CM (1998). Modeling of storm responses in conduit flow aquifers with reservoirs, J. Hydrol., 208, 82-91. Harbaugh AW (2005) MODFLOW-2005, the U.S. Geological Survey modular groundwater model— The ground- water flow process: U.S. Geological Survey Techniques and Methods 6-A16. Harbaugh AW, Banta ER, Hill MC & McDonald MG (2000). MODFLOW-2000, The U.S. Geological Sur- vey modular ground-water model—User’s guide to modularization concepts and the ground- water flow process: U.S. Geological Survey Open-File Report 00–92, 121 p. Hartmann A, Goldscheider N, Wagener T, Lange J & Weiler M (2014). Karst water resources in a changing world: review of hydrological modeling approaches. Rev Geophys 52:218–242. Hosseiny Sohi SM, Koch M & Ashjari J (2014). Construction Effects of the Karun 4 Dam, Iran, on the Groundwater in the Adjacent Karstic Aquifer. ICHE 2014, Hamburg - Lehfeldt & Kopmann (eds). Hosseiny Sohi SM, Koch M & Ashjari J (2017a). Evaluating permeability and groutability at Karun 4 dam Iran using Lugeon values and grout Take. Symposium Proceeding of 85th Annual Meeting of International Commission on Large Dams, Prague, July 2017. Hosseiny Sohi SM, Koch M & Ashjari J (2017b). Monitoring of the seepage controlling system of the Karun 4 dam in Iran. Symposium Proceeding of 85th Annual Meeting of International Commission on Large Dams, Prague, July 2017. Houlsby AC (1990). Construction and design of cement grouting. A guide to grouting in rock foundations. New York, Wiley (Wiley series of practical construction guides). Huntoon P W (1995). Is it appropriate to apply porous media groundwater circulation models to karstic aquifers El-Kadi, A. I., ed., Groundwater models for resources analysis and management, 1994 Pacific Northwest/Oceania Conference, Honolulu, HI, p. 339-358. Jaeger Ch (1979). Rock mechanics and engineering. Cambridge University Press. Jeannin P-Y (2001). Modeling flow in phreatic and epiphreatic karst conduits in the Holloch Cave (Muotatal, Switzerland). Water Resour. Res., 37, 191-200. Keeler RR & Zhang YK (1997). Modeling of groundwater flow in a fractured-karst aquifer in the Big Springs Basin, Iowa. GSA Abs. with Programs, 29, 4, 25. 90 Kiraly l (1998). Modeling karst aquifers by the combined discrete channel and continuum approach. Bulletin du Centre d’hydrogeologie, Neuchatel, 16, 77-98. Knill, JL (1972). Assessment of reservoir feasibility. Quart. J. Eng. Geol, 4, 355-372. Kuniansky EL & Shoemaker WB (2009). Conduit Flow Process (CFP) for MODFLOW- 2005: in Proceedings 15th International Congress of Speleology, Kerrvile, Texas, USA, July 19- 26, 2009; William B. White, ed.: Greyhound Press, Huntsville, Alabama: pp. 1568-1574. Mahab Ghods Consulting Engineering Co. (2002). Final engineering geological report (phase 2), Karun 4 dam project. Mahab Ghods Consulting Engineering Co. (2010). Final excavation and injection report of the Karun 4 dam grout curtain. Mahab Ghods Consulting Engineering Co. (2015). Monitoring report of Karun 4 dam. McDonald MG & Harbaugh AW (1988). A modular three-dimensional finite-difference ground-water flow model: U.S. Geological Survey Techniques of Water Resources Investigations, book 6, chap. A1, 548 p. Milanović PT (2004). Water resources engineering in karst. CRC Press, Boca Raton. Mohrlok U & Teutsch G (1997). Double continuum porous equivalent (DCPE) versus discrete modeling in karst terranes. In: Gunay G, Johnson AI (eds) Karst waters and environmental impacts. In: “Proceedings 5th Int Symp Field Sem Karst Waters and Environmental Impacts, Turkey. Balkema, Rotterdam, pp. 319–326. Neuman SP (2005). Trends, prospects and challenges in quantifying flow and transport through fractured rocks. Hydrogeology Journal 13 (1), 124 - 147. Quinlan JF & Ewers RO (1985). Groundwater flow in limestone terrains: Strategy rationale and procedure for reliable and efficient monitoring of groundwater quality in karst areas. In: Proceed. 5th national symposium and exposition on Aquifer Restoration and Groundwater Monitoring, The Fawcett Center, Columbus, Ohio. pp 197-234. Reimann T, Geyer T, Shoemaker WB, Liedl R & Sauter M (2011). Effects of dynamically variable saturation and matrix-conduit coupling of flow in karst aquifers, Water Resources Research, v. 47, W11503, doi:10.1029/2011WR010446. Reimann T & Hill ME (2009). MODFLOW-CFP: A New Conduit Flow Process for MODFLOW-2005: Ground Water, v. 47, p. 321-325. Sadd Tunnel Pars Company (2015). Evaluation of Grouting Curtain of Karun 4 dam. Sahuquillo A (1985). Spanish experience in karst water resources. In: Proceedings, international symposium on karst water resources, Ankara, Turkey, IAHS Publication 161, pp 133–147. Saller SP, Ronayne MJ & Long AJ (2013). Comparison of a karst groundwater model with and without discrete conduit flow: Hydrogeology Journal, v. 21, no. , p. 1555-1566. Sauter M (1992). Quantification and forecasting of regional groundwater flow and transport in a karst aquifer (Gallusquelle,Malm,SW Germany): Tübinger Geowissenschaftliche Arbeiten, Ser.C, 13, 150p. 91 Shoemaker W B, Kuniansky EL, Birk S, Bauer S & Swain ED (2008). Documentation of a conduit flow process (CFP) for MODFLOW-2005, U.S. Geol. Surv. Tech. Methods, Book 6, Chap. A24. Uromeihy A (2000). The Lar Dam; an example of infrastructural development in a geologically active karstic region. Journal of Asian Earth Sciences, Volume 18, Issue 1, pp 25–31. Vennard JK & Street RL (1975). Elementary fluid mechanics (5th ed.): New York, John Wiley, 740 p. White WB (1999). Conceptual models for karstic aquifers. In: Karst Modeling, Karst Waters Inst. Spec. Publ., Vol. 5, A. N. Palmer, M. V. Palmer and I. D. Sasowsky (eds), pp. 1-16, Karst Waters Inst., Charles Town, W. Va. Winston RB (2009). ModelMuse: a graphical user interface for MODFLOW-2005 and PHAS T. US Geological Survey, Reston, VA. Worthington SRH, Ford DC & Beddows PA (2000). Porosity and permeability enhancement in unconfined carbonate aquifers as a result of solution, In “Speleogenesis: Evolution of Karst Aquifers” A. B. Klimchouk et al. (eds)., pp. 463-472, Natl. Speleology Soc., Huntsville, Al. Yurtsever Y & Payne R.B (1986). Mathematical models based on compartmental simulation approach for quantitative interpretation of tracer data in hydrological systems, Proc. 5th Intl. Symposium of Underground Water Tracing, IGME, Athens, Greece, pp. 341-353. Zhang YK, Bai EW, Libra R, Rowden R & Liu H (1996). Simulation of spring discharge from a limestone aquifer in Iowa, Hydrogeol. J.., 4, 41-54. 92 Chapter 5 - Investigation of different Scenarios of Water Leakage through Foundation and Abutments of Karun 4 Dam Abstract Despite all recent achievements in Karst aquifers’ research, the numerical modelling of groundwater flow processes in such aquifers has remained complicated and problematic, since the precise determination of the conduits network within a particular karst aquifer requires a significant amount of local data, observations by site experts, and sound application findings from laboratories/academic institutes researchers. The objective of the present study is to investigate different scenarios of water leakage through the karstic dam foundation and abutments of Karun 4 dam in Iran. To this avail, classical groundwater flow in the porous matrix medium and laminar or turbulent pipe flow in the discrete conduits are simulated with the recently developed numerical groundwater flow model MODFLOW-CFP. A total of 97 scenarios, divided into three categories, are simulated and the computed amount of downstream water escape compared with that of the calibrated transient reference model scenario (RS): In the first category, for each of the 15 originally recognized karstic faults, represented by calibrated pipes in the reference CFP-model, four different diameters: 20, 50, 200 and 500 % of the corresponding reference diameter of the pipe, are analyzed. In the second category, 10 new hypothetical pipes are embedded at the location of extra faults with three scenarios for each pipe. In the third category several adjusted pair-pipes that represent continuous faults on the up- and downstream of the grout curtain are connected to each other. Results of the first group of scenarios identified the most critical fault, namely, F4-Up, whose pipe characteristics have the largest effect on the downstream water escape. Results of the second scenario category emphasize that some potential faults at the border of the adjacent As1-g and As1-i karstic formation units may lead to an additional water discharge of nearly 1% relative to that of the RS. Finally, the connected-pipes scenarios results demonstrate the effectiveness of the grout curtain in impeding downstream water leakage even if potential karstic faults close to the critical fault F4-Up may open up further. Key words: Carbonate aquifer and karst, Conduit Flow Process (CFP), Numerical modelling, Water leakage scenarios, Karun 4 dam, Iran 93 5.1 Introduction Carbonate aquifers with karst potential and ensuing large fissures/conduits supply a huge amount of drinkable water in terms of groundwater of our blue globe (Ford and Williams, 2007). On the other hand, because of their complexity and the inability to describe such aquifers conceptually as a Darcy porous medium, as well as a general lack of reliable data, fractured carbonate aquifers are considered to be one of the most challenging research topics in the general field of hydrogeology, last but not to the least (Hartmann et al., 2014). Furthermore, when it comes to the construction of a high concrete dam in an appropriate narrow valley with its abutments constructed over carbonate limestone, additional concerns arise, because increasing water level adjacent the carbonate aquifer could cause additional karstification, which may lead to proceeding development of fissures and fractures, that may expand to conduits and even caverns and caves (Milanović, 2004; Ford and Williams, 2007). The final result of this process could be water leakage through the dam foundations and abutments which, consequently, might prevent a complete filling of the dam reservoir.. This is exactly what happened at the Camaras dam in Spain (Jaeger, 1979) and the Lar dam in Iran (Uromiyehi, 2000). The classification of Karst aquifers based on their hydraulic specifications has been presented by Quinlan and Ewers (1985) and they are as follows: - Mature system with conduits that have a high hydraulic conductivity, and which has been formed by fracturing and faulting processes or by the dissolution of conduits that contains essentially turbulent pipe flow as the pore openings are more than 1 cm wide. - Dispersion system which has a conduit matrix, with openings less than 1cm wide and where groundwater tends to behave as laminar flow. - Hybrid system which includes both conduit and matrix systems. Most karst aquifers can be classified as a hybrid system. Hence, to model karst aquifers, the special porosity structure of such a media, namely, the double continuum porosity must be considered. Barenblatt et al. (1960) were the first to present this duality concept, which was investigated further and its applicability to Karstic aquifers demonstrated by Atkinson (1977), Gunn (1985), Yurtsever (1986), Huntoon (1995), Zhang (1996), Keeler (1997), Kiraly (1998), 94 Halihan & Wicks (1998), White (1999), Jeannin (2001), Neumann (2005), and Ghasemizadeh et al. (2012). Repudiated numerical models to simulate groundwater flow in a karst media are the following: - Double Continuum Porous Equivalent model (Sauter,1992), - CAVE (Carbonate Aquifer Void Evolution) model (Clements et al.,1996), - GW (groundwater) model (Cornaton and Perroche, 2002), and, - CFP (Conduit Flow Process) for MODFLOW-2005 (Shoemaker et al.2008). The differences between these models are their methodology to define the conduit medium. In the Double Continuum Porous Equivalent model, the conduit medium is presented by an equivalent porous medium where the classical Darcy law is applicable. In the CAVE and CFP models the conduits are adjusted by a pipe network in which turbulent flow occurs. The GW- model is essentially the same as the CAVE and CFP models, but uses a finite element- instead of finite difference method. “As MODFLOW is so widely circulated and used in hydrogeology world, efforts are being made to adapt the code for Karst conduits” (Ford and Williams, 2007). One year after this statement, the MODFLOW-CFP module as an add-on to the MOFLOW-2005 code (Harbaugh, 2005).was released by Shoemaker and colleagues of the U.S. Geology Survey. Many researchers applied CFP since that time and evaluated its validity in comparison to either previous similar codes or laboratory tests. Hill et al. (2008), for instance, employed CFP to investigate the specifications of a karst aquifer underlying west-central Florida. Reimann et Hill (2009) and Kuniansky and Shoemaker (2009) described the characteristics of CFP. Ashok and Sophocleous (2009) evaluated the recent MODFLOW-developments for more advanced groundwater flow modelling. Reimann et al. (2011a) inspected the effects of variable saturation and matrix-conduit coupling of flow in karst aquifers, and Reimann et al. (2011b) subsequently published a modification to CFP. Saller et al. (2013) compared a karst groundwater model flow with and without discrete conduit flow. Gallegos et al. (2013) applied CFP to simulate flow in Karst aquifers at laboratory- and sub-regional scales. Assari and Mohammadi (2017) used CFPM2 (Mode 2 of CFP) and stochastic simulations to assess flow paths in the Karst aquifer of the Asmari formation in southwest Iran . 95 The subject of the present study is the Karun 4 dam in Iran which is constructed mainly on a limestone carbonate aquifer, namely, the Asmari 1 formation. Due to this porous/conduit limestone specifications, the risk of water leakage through the faults and fractures should be considered seriously. Although to deal with this problem, a grout curtain was drilled and injected through grouting galleries in the dam foundations and abutments and ground water levels, monitored by a set of piezometers and observation wells, the need of comprehensive investigation via a numerical model is necessary. This study relies on the results of the recent works of the author(s) regarding Karun 4 dam groundwater flow. The chemical investigations of water in dam reservoir and downstream indicated a significant arising of the ground water level in the area after dam impounding (Hosseiny Sohi et al., 2014). Using the method proposed by Ewert (1985), the permeability of the grout curtain was estimated by means of Lugeon tests- and cement-take results at the dam abutments and foundation (Hosseiny Sohi et al., 2017a). Afterwards, changes in groundwater heads were monitored and analyzed in standpipe piezometers of the grouting galleries and observation wells located in different places around the dam area (Hosseiny Sohi et al., 2017b). Finally, a transient numerical groundwater flow model using the CFP- module for MODFLOW- 2005 was designed and calibrated to simulate the groundwater flow in the dual porosity media of the Karun 4 dam foundations and abutments (Hosseiny Sohi et al., 2018). The CFP- MODFLOW- model can simulate successfully laminar flow in the matrix media and laminar/turbulent flow in the conduits and will be used in the present analysis to study the different scenarios regarding water leakage as follows: - Adjusting the pipe diameters in the model and investigating its impact on the water escape outputs, - Inserting and adjusting hypothetical pipes in the model to find out new possible water leakage paths, - Connecting the adjusted pair-pipes to each other, to represent the continuous faults on the up- and downstream sections of the grout curtain, to study the possibility of conduits development and its effect on the efficiency of the grout curtain at the corresponding places. The results of these scenario-simulations should provide hints to identify the concerning and critical water leakage paths and which need further in-situ monitoring and observations. 96 The present numerical MODFLOW/CFP- model is a novel approach in this topic of engineering hydrogeology, as it couples laminar Darcy- flow in the porous matrix medium with pipe turbulent flow and applies this methodology to a very specific dam site where changes of the groundwater level of more than 150 meters occur in a very short distance up- and downstream of the reservoir. 5.2 Study area 5.2.1 Site specifications and geological setting The Karun 4 dam is located at geographical latitude 50° 28´ 19´´and longitude 31° 35´ 50´´, close to the source of the Karun River in Chaharmahal and Bakhtiari province, southwest Iran, upstream of Khuzestan province, where other important national dams on the Karun river have been constructed and are under operation. The dam reservoir with an area of 29 km2 connects the Armand and Bazoft rivers in the upstream of Karun river. The total reservoir volume reaches 2.2 billion cubic meters for the maximum normal water level of 1028 m.a.s.l. The dam’s purpose is to regulate the Karun river’s discharge and to control the frequent destructive floods of the Karun River, but also to generate about 1GW hydroelectric energy (Figure 5.1). The Karun 4 dam is in the tectonically folded belt of the Zagros Mountains within a canyon that is located in the southwestern flank of the Kuh-Sefid anticline which has a NW-SE trend. Most of the reservoir is situated on Cretaceous to Miocene limestone and marly limestone (Haghnejad et al., 2014). The geological formations around the dam site consist of carbonate layers of the Asmari formation (As) which can be divided further into lower (As1), middle (As2) and upper (As3) units and marlstone and marly limestone of the impervious Pabdeh formation (Pd) upstream of the dam site (see Figure 5.1). The geology of the dam site is affected by many faults. The lower unit of the Asmari formation (AS1) consists of limestone, marly limestone and marlstone with some porous limestone (Hosseiny Sohi et al, 2014), and is subdivided from upstream to downstream into 10 more detailed subunits AS1-a to AS1-j, as described in more details in Table 5.1. 97 Figure 5.1. Plan view of the Karun 4 dam in Iran Table 5.1. Lithological specifications of Asmari 1 subunits Name Description AS1-a Limestone (about 90%) gray to brownish gray, thick to very thick, strength to very strength, with no considerable karstic phenomena. There is some surficial erosion which is related to marly limestone (about10%) AS1-b Limy marlstone mostly covered by colluviums and slope wash. This rock is gray to dark green in natural color, weak to relatively strength and with no karstic features AS1-c Limestone, thick to very thick bed, gray to light brown with no particular karstic phenomena. AS1-d More or less like AS1-b AS1-e Limestone with some (in unweathered locations) gray to blue outcrops (marly limestone) very thick layers and without considerable karstic phenomena AS1-f This unit commences at As1-e and continues to the first porous key bed. Composed mostly with limestone having slight to medium karstic features (vug) and 10% porous limestones. AS1-g Porous limestone which are solution vugs gray to light brownish gray in natural colour. At the right abutment the fresh outcrop is gray. It is the key porous bed. AS1-h Limestone with some grey to blue outcrops (probably marly limestone very thickbed. There are more or less karstic vugs at the top of this sub-unit AS1-i This sub-unit is more or less similar to AS1-g, so it is considered as key bed. AS1-j Lime stone with marly limestone. the Marl limestone is frequent in the lower part, thick to very thick, no important karstic features is seen in this sub-unit 98 Figure 5.2. Grout curtain expanded cross section 5.2.2 Seepage controlling system To seal the dam foundations and abutments, a very large grout curtain was designed and constructed. The construction was done through the excavation of 5 series of galleries in both left and right abutments at different elevations from the dam crest elevation to the bottom and, finally, the foundation gallery, called the 806-gallery (see Figure 5.2). From inside these galleries, the grouting boreholes are drilled, examined and injected to have an integrated virtual wall or barrier, called grout curtain against water leakage (Houlsby, 1990). Some of the primary drilled holes were selected as exploratory boreholes and the rock permeability tests and geological logs were carried out within these. They were then considered in the grouting phase as normal drilled boreholes and the amounts of their cement Take were reported (Hosseiny et al, 2017b). Based on the results obtained from these exploratory boreholes, the grout curtain in two lines was designed and constructed to sew/stitch the dam structure to the Pabdeh impermeable formation in the upstream section of the dam. The seepage controlling system comprises the grout curtain and the drainage curtain, wherefore the latter consists of boreholes drilled at 40 meters’ horizontal distance toward downstream in two galleries in the left and right abutments and allows not only to detect and to explore the still leaking water through the grout curtain, but also aids to reduce the uplift pressure and so prevents the breaking of mass rock downstream. To monitor the dam’s hydraulic behavior, a network of instruments is installed in the dam body, grouting galleries and at different stations of the dam site. The most important 99 instruments of the seepage controlling system are the standpipe piezometers and observation wells. The former, also known as Casagrande piezometers, are installed in the grouting galleries upstream and downstream of the grout curtain, whereas the observation wells are drilled in the different points of the dam body, abutments, the drainage curtain and downstream of the dam site. 5.3 Mathematical and numerical methodology 5.3.1 Theoretical basics of the MODFLOW-CPF dual porosity model The numerical groundwater flow model to simulate flow process in a karstic dual porosity media with discrete faults/fissures is the Conduit Flow Process (CFP) (Shoemaker et al., 2008) module for MODFLOW-2005 (Harbaugh, 2005). “CFP has the ability to simulate turbulent or laminar groundwater flow conditions by coupling the traditional groundwater flow equation with formulations for a 1-dimensional discrete network of cylindrical pipes (Mode 1, CFPM1).” (Shoemaker et al., 2008) The governing equation applied to the porous medium is the general groundwater equation derived from Darcy law (1856) and the mass-conservation law. Darcy’s law for a one-dimensional flow in a prism of porous material (Figure 5.3) is defined as (McDonald and Harbaugh, 1988): 𝑄 = 𝐾𝐴(ℎ1−ℎ2) 𝐿 (5.1) where K [LT-1] is the hydraulic conductivity of the matrix in the direction of flow A [L2] is the cross-sectional area perpendicular to the flow ℎ1 − ℎ2 [L] is the head difference across the prism parallel to flow L [L] is the length of the prism parallel to the flow path. From the continuity equation could be deduced that the sum of all flows into and out of the cell must be equal to the rate of change in storage within the cell. Under the assumption that the density of ground water is constant, the continuity equation expressing the balance of flow for a cell is (McDonald and Harbaugh, 1988): ∑ 𝑄𝑖 = 𝑆𝑠 𝜕ℎ 𝜕𝑡 ∆𝑉 (5.2) 100 Figure 5.3 Darcy’s law for a prism of porous material (from McDonald and Harbaugh, 1988) where Qi [L3T-1] is the flow rate across face i of the cell (= positive for inflow and negative for outflow) Ss [L-1] is specific storage, ΔV [L3] is the volume of the cell, and ℎ [L] is the potentiometric head; 𝑡 [T] is time. The term on the right-hand side is equivalent to the net change of volume of water taken into/from the storage over a time interval 𝜕t given a change in head of 𝜕h. Combining equations (5.1) and (5.2) together, the general governing differential equation can be derived (Anderson et al., 2015), which represent three-dimensional transient groundwater flow for heterogeneous and anisotropic conditions (McDonald and Harbaugh, 1988): 𝜕 𝜕𝑥 (𝐾𝑥𝑥 𝜕ℎ 𝜕𝑥 ) + 𝜕 𝜕𝑦 (𝐾𝑦𝑦 𝜕ℎ 𝜕𝑦 ) + 𝜕 𝜕𝑧 (𝐾𝑧𝑧 𝜕ℎ 𝜕𝑧 ) ± 𝑊 = 𝑆𝑠 ( 𝜕ℎ 𝜕𝑡 ) (5.3) where 𝐾𝑥𝑥, 𝐾𝑦𝑦, 𝐾𝑧𝑧 [LT -1] are the values of hydraulic conductivity along the x, y, and z axes, respectively; 𝑊 [T-1] is a volumetric source flux per unit volume; 𝑆𝑠 [L -1] is specific storage. This equation is applied for the numerical modelling by MODFLOW packages and assumes that ground water is incompressible with constant density and viscosity. The equation describing flow in the conduits is the Darcy-Weisbach equation for both laminar and turbulent flow (Weisbach,1845; Brown, 2003): 101 ∆ℎ = ℎ𝐿 = 𝑓 ∆𝑙 𝑑 𝑉2 2𝑔 (5.4) where: ∆ℎ or ℎ𝐿 [L] is the head loss measured along the pipe of length ∆l [L], 𝑓 [unitless] is the friction factor, 𝑑 [L] is pipe diameter, 𝑉 [LT-1] is the mean velocity, and 𝑔 [LT-2] is the gravitational acceleration constant. The exchange discharge between the matrix medium described by MODFLOW-2005 and the conduit is computed through the following equation in CFP (Shoemaker et al., 2008): 𝑄𝑒𝑥 = 𝛼𝑗,𝑖,𝑘(ℎ𝑖𝑛 − ℎ𝑗,𝑖,𝑘) (5.5) where 𝑄𝑒𝑥 [L3T-1] is the volumetric exchange flow rate, 𝛼𝑗,𝑖,𝑘 [L 2T-1] is the pipe conductance at MODFLOW cell j,i,k ℎ𝑖𝑛 [L] is the head at the pipe node n located at the center of the MODFLOW cell, ℎ𝑗,𝑖,𝑘 [L] is the head in the encompassing MODFLOW cell j,i,k, and 𝛼𝑗,𝑖,𝑘 is computed by CFP as: 𝛼𝑗,𝑖,𝑘 = ∑ 𝐾𝑗,𝑖,𝑘𝜋𝑑𝑖𝑝 1 2 (Δ𝑙𝑖𝑝𝜏𝑖𝑝) 𝑛𝑝 𝑖𝑝=1 𝑟𝑖𝑝 (5.6) where 𝐾𝑗,𝑖,𝑘 [LT -1] is the conduit wall permeability (conductivity) of MODFLOW cell j,i,k, 𝑑𝑖𝑝 [L] is the diameter of pipe ip connected to node in, 𝛥𝑙𝑖𝑝 [L] is the straight-line length between nodes of pipe ip connected to node in, 𝜏𝑖𝑝 [unitless] is the tortuosity of pipe ip connected to node in, 𝑟𝑖𝑝 [L] is the radius of the pipe. 102 Figure 5.4 Boundary conditions layout of the numerical model As the developers pointed out, CFPM1 can represent dissolution or biological burrowing features in carbonate aquifers and voids in fractured rock. Thus, Mode 1 is employed here to simulate the karstified limestone formation of Asmari 1 in the Karun 4 dam area. 5.3.2 Set-up of the CPF- model 5.3.2.1 General layout and boundary conditions Using various data, such as topography and geological maps, piezometers and observation wells measurements, reservoir water level, drainage data and layout of dam body, grout curtain and drainage curtain, mostly received from Iran Water and Power development Company (IWPC), ModelMuse, a graphical user interface of the U.S. Geological Survey (USGS) (Winston, 2009) is applied to set up the model. The specifications of the model are briefly as follows: The grid cell size is 5x5 m2, so that the whole study area includes 170 rows and 150 columns in north-south and west-east directions, respectively, i.e. a total of 25500 cells in each horizontal plan view, however where most of them are not active cells. The model includes four layers, with their elevations, types and active cells listed in Table 5.2 . The boundary conditions used are displayed in Figure 5.4. 103 Table 5.2. Layer specifications assigned to the aquifer model Name Elevation (m.a.s.l) Type of Aquifer Number of active cells Bottom Top Layer 1 (top layer) 974 1032/topography unconfined 6,246 Layer 2 890 974 confined 7,747 Layer 3 848 890 confined 9,248 Layer 4 (bottom layer) 700/Pabdeh 848 Confined 10,217 Hydraulic conductivity (K) variations of each layer are derived from different zones of Asmari1 subunits, Pabdeh- and Asmari2 formations. Figure 5.5 shows the corresponding zonation of K for layer 2 as a typical map for assigning the K-value zones in the model. Later, in the transient state modeling, the same zones are assigned to storage specifications. Dam body, grout curtain and drainage curtain layout are imported as shape files from ArcGIS into ModelMuse. The observed data of piezometers and observation wells (see Figure 5.6) are used in the Head Observation Package (HOB) of MODFLOW-2005 to compute the residuals between observed and simulated ground water heads. Figure 5.5 Geological map at layer 2 including Asmari 1 subunits 104 Figure 5.6 Layout of piezometers and observation wells The model is calibrated afterwards for steady- and transient state. In transient state, time series data of four years and nine months, from June 2012 up to March 2017 (=1727 days), are divided into 38 stress periods for the calibration phase and 21 stress periods for the validation phase. Moreover, each stress period is made up in the present model of one time step, with a slightly variable length of around 30 days, on average. 5.3.2.2 Incorporation of karst conduits and final calibration of the CPF-model After adjusting the different calibration parameters of the groundwater model for the classical Darcy porous medium (see Sections 4.4.1 and 4.5.1), the pipes are added to the CFP model to represent the conduits (turbulent or laminar flow). The layout of the possible pipes was obtained from the exiting faults in the geological maps. The corresponding CFP parameters to be assigned to each pipe are the following: Pipe diameter, Tortuosity, Roughness height, Lower- and Higher critical Reynolds number, Conduit wall permeability (hydraulic conductivity) and Elevation of pipe. However, the calibration parameters of the pipes are only its diameter and conduit wall permeability. Other possible parameters have no significant impact on the calibration process. Hence, the following inherent values are assigned to all pipes: 105 Tortuosity: 1.0, Roughness height: 0.01, Lower critical Reynolds number: 2, Higher critical Reynolds number: 10 (Shoemaker et al., 2008), Elevation of pipe: 895, 850 and 800 m.a.s.l to layer 2, 3 and 4 respectively. With the above-mentioned assumptions, trial- and error calibration is done for each pipe to find out the most appropriate diameter and permeability K. In both pipe diameter D- and K- calibration, the only criteria to decide to keep the corresponding value, or even the fault as a pipe, was its effect on the rmsr of the model, i.e. reducing it. Therefore, as not all candidate faults retrieved from the geological maps could be calibrated as discussed above, they were discarded from the CFP- model. Consequently, the finally defined 15 pipes of the model illustrated in Figure 5.7 describe only those faults whose characteristics can be most likely expressed as a discrete conduit within the matrix medium in a karstified formation to be simulated with the CFPM1- model mode (see Section 4.4.3). The finally calibrated parameters, namely, diameters and K-values of the pipes and the fault they represent, are listed in Table 5.3. Table 5.3. Finally calibrated pipe specifications (K- and D- values) of the CFP- model and association of the various pipes with the real existing faults. Elevation (m.a.s.l) Pipe K (m/day) Diameter (m) Length (m) Start Node Final Node Fault Location 848 (layer 4) 1 0.005 0.02 48 1 11 F4-Up Upstream 2 0.005 0.1 211 12 71 F7 Downstream 3 0.00005 0.01 116 72 91 F13 Upstream 4 0.09 0.005 106 92 113 F13 Downstream 5 0.1 0.02 310 114 184 F4-Up Downstream 890 (layer 3) 6 0.1 0.1 93 185 207 F2 Upstream 7 0.12 0.079 138 208 246 F3 Downstream 8 0.01 0.007 47 247 257 F4-Up Upstream 9 0.05 0.01 299 258 325 F4-Up Downstream 10 0.09 0.05 283 326 389 F4-Down Downstream 932 (layer 2) 11 0.0005 0.01 49 390 401 F2 Downstream 12 0.00001 0.05 43 402 413 F2 Upstream 13 0.00007 0.01 91 414 440 F3 Upstream 14 0.12 0.1 171 441 488 F15 Upstream 15 0.00001 0.05 101 489 517 F14 Upstream The criterion of succeeding a good calibration and validation is the reduction (whatever is possible) of the root mean square residual (rmsr) between the observed and simulated ground water heads in the piezometers and observation wells. The results are listed in Table 5.4 for both steady- and transient state, wherefore the latter includes calibration and validation phases. 106 Figure 5.7 Layout of the pipes at a) layer 2, b) layer 3, c) layer 4, derived from the location of the detected faults from geological maps, d) faults plan view 107 Table 5.4. Statistical results of the calibration the model for steady-state and transient (calibration and validation) states State Rmsr(m) R2 Steady 1.44 0.998 Transient (calibration) 2.28 0.98 Transient (validation) 2.60 0.97 Figures 5.8 and 5.9 show the transient state results of groundwater heads in calibration and validation phases. Figure 5.8 Transient state simulated- versus observed heads for a) calibration- and b) validation phase Figure 5.9 Transient state residuals in calibration- (left) and validation phase (right) 108 5.4 Simulation scenarios of pipe adjustments for the representation of faults 5.4.1 Overview of scenario specifications and definition of the reference scenario The present calibrated model whose specifications is described in the previous Section 5.3 enables one to investigate different scenarios of water leakage within the karstified carbonate aquifer of Asmari 1 on the foundation and abutments of the Karun 4 dam. In this section, 97 scenarios separated in three categories as following are analyzed and compared with the reference scenario (RS) (see below): - Changing of adjusted pipes diameters in the model to study their impact on the water escape outputs (60 scenarios), - Embedding of hypothetical pipes in new locations in the model area to find out new possible water leakage paths through karst conduits/faults (28 scenarios), - Connecting adjusted pair-pipes to each other to represent continuous faults on the up- and downstream section of the grout curtain and so to study the possibility of new conduits’ developments and their effects on the efficiency of the grout curtain at the corresponding places (9 scenarios). The model is run for all these scenarios and the results of the different parameters, like water discharge of pipes, drainage, water discharge through the downstream boundaries (constant head OUT), storage in the aquifer are extracted and compared with the results of the Reference Scenario (RS), namely, the transient state model (see Table 5.5). Table 5.5. Cumulative volume budget in (m3) for the entire model at the end of the simulation period (of 1727 days) for the Reference Scenario (RS) IN OUT Storage 28,250,671 Storage 10,649,783 Constant Head 151,494,614 Constant Head 90,914,165 Drains 0 Drains 78,182,870 Pipes 1,482,481 Pipes 1,482,481 Total IN 181,227,767 Total OUT 181,229,300 The simulation results of these three scenario-categories listed above are discussed in the subsequent sub-sections. 109 5.4.2 Scenarios with changing pipe diameters To study the impact of the pipe diameter variation on the whole model outcome, different scenarios classified in the following order are simulated: For each of the 15 original pipes, four different diameters: 20, 50, 200 and 500 % of the corresponding reference diameter of a pipe, i.e. a total of 60 diameter scenarios, are analyzed. 5.4.2.1 Pipe discharge results Table 5.6 demonstrates the percentile changes in the pipe’s water discharge (see column: Pipes). The differences for the 11 scenarios are more than 10%, among them 6 scenarios gain less volume of water in the pipes (negative values), and for 5 scenarios more volume of water flows through the pipes (positive values). Meanwhile, when increasing the pipe diameter to 500%, the model does not converge for pipe 7 at elevation 890m. The scenarios with the highest percentile changes of pipes’ water discharge are presented in Table 5.7 and as a bar diagram in Figure 5.10. Obviously, augmenting the pipe diameter D will increase (positive sign) the pipe’s discharge volume, while reducing it will decrease (negative sign) it. Increasing discharge amounts are occurring for pipes 9 and 10 in layer 3 and pipe 5 in layer 4. All of them are located at the left side and downstream of the dam and grouting curtain, and they are representatives of faults F4 Up (pipe 5 and pipe 9) and F4 Down (pipe 10) (Mahab Ghods report, 2002). In pipe 5, increasing the pipe diameter D by a factor ranging from two to five increases the pipe’s volume discharge from 24% to 63 %, whereas in pipe 10 located close to the former, but in a layer above, the increase varies only from 22% to 25 %, due to the fact that D of the first pipe is 0.02m and that of latter 0.05 m. The other six pipes experience a decline of their volume discharge of around 20 % which means that they will be sealed to some extent. 110 Table 5.6. Changes (%) for the different pipe diameter scenarios relative to the reference scenario RS Number Hydraulic conductivitiy (m/day) Diameter (m) Diameter % Diameter (m) rmsr STORAGE CONSTANT HEAD DRAINS PIPES 2.43880 17600888.7 90914165.41 78182870 1482481.402 1 20% 0.004 2.43880 0.00059 -0.00038 -0.00023 -0.30526 2 50% 0.01 2.43879 0.00142 0.00003 -0.00005 -0.11503 3 200% 0.04 2.43880 0.00046 -0.00006 -0.00002 0.00517 4 500% 0.1 2.43880 0.00058 -0.00006 -0.00001 0.00491 5 20% 0.02 2.43960 0.00051 0.00242 0.00327 -1.89644 6 50% 0.05 2.43883 -0.00026 0.00009 0.00010 -0.04635 7 200% 0.2 2.43880 0.00072 0.00000 -0.00001 0.00261 8 500% 0.5 2.43882 0.00007 -0.00001 0.00000 0.00247 9 20% 0.002 2.45681 0.30090 0.09343 0.01161 0.14798 10 50% 0.005 2.43881 0.00131 -0.00011 -0.00006 -0.00702 11 200% 0.02 2.43874 -0.00152 0.00340 0.00092 -0.03188 12 500% 0.05 2.43881 -0.00038 -0.00012 -0.00002 0.00042 13 20% 0.001 2.43881 -0.00044 -0.00016 -0.00019 -0.05816 14 50% 0.0025 2.43877 0.00297 -0.00039 -0.00018 -0.05827 15 200% 0.01 2.43885 0.00008 0.00044 0.00139 0.41710 16 500% 0.025 2.43942 -0.00101 0.00278 0.00949 3.80870 17 20% 0.004 2.44237 -0.00347 0.00516 -0.02248 -5.89081 18 50% 0.01 2.44201 -0.00242 0.00482 -0.01930 -4.31111 19 200% 0.04 2.42644 0.01055 -0.03044 0.09431 23.62304 20 500% 0.1 2.43435 0.15172 0.00945 0.14833 63.46099 21 20% 0.02 2.45485 0.27912 0.05122 0.00255 -19.12184 22 50% 0.05 2.43848 -0.00967 -0.00973 -0.00212 -4.24536 23 200% 0.2 2.43882 -0.00160 0.00018 0.00004 0.13540 24 500% 0.5 2.43882 -0.00041 0.00014 0.00004 0.13975 25 20% 0.0158 2.47565 -0.01088 -0.03368 0.01441 -20.80702 26 50% 0.0395 2.45537 -0.00633 -0.02268 0.00924 -12.78189 27 200% 0.158 2.43739 -0.00603 0.00272 -0.00135 1.60174 28 500% 0.395 29 20% 0.0014 2.43879 0.00180 -0.00009 -0.00003 -0.08927 30 50% 0.0035 2.43881 -0.00115 -0.00022 -0.00005 -0.08380 31 200% 0.014 2.43883 0.00080 0.00024 0.00006 0.27013 32 500% 0.035 2.43883 0.00151 0.00035 0.00009 0.35193 33 20% 0.002 2.44466 0.04335 0.14216 0.01628 0.45120 34 50% 0.005 2.44367 0.04352 0.13178 0.01468 0.63207 35 200% 0.02 2.45549 0.29062 0.09337 0.01184 1.40776 36 500% 0.05 2.44480 0.06696 0.17313 0.00392 14.90144 37 20% 0.01 2.46260 -0.03123 -0.04825 0.00733 -18.30766 38 50% 0.025 2.45719 -0.01920 -0.04149 0.00790 -14.97392 39 200% 0.1 2.43360 0.03810 0.08625 -0.02541 22.48133 40 500% 0.25 2.43455 0.04818 0.09880 -0.02913 25.33600 41 20% 0.002 2.43880 0.00011 -0.00006 -0.00002 0.00048 42 50% 0.005 2.43881 -0.00078 -0.00005 -0.00001 -0.00811 43 200% 0.02 2.43881 0.00059 -0.00004 -0.00001 0.00147 44 500% 0.05 2.43881 0.00127 -0.00011 -0.00003 0.00402 45 20% 0.01 2.43881 0.00109 -0.00013 -0.00003 0.00443 46 50% 0.025 2.43881 0.00004 -0.00004 -0.00001 0.00053 47 200% 0.1 2.43882 0.00114 -0.00011 -0.00003 0.00135 48 500% 0.25 2.43880 0.00031 -0.00001 -0.00001 0.00235 49 20% 0.002 2.45628 0.28360 0.09056 0.01142 0.16028 50 50% 0.005 2.43881 0.00066 -0.00007 -0.00003 -0.01132 51 200% 0.02 2.43881 -0.00035 -0.00008 -0.00002 0.00045 52 500% 0.05 2.43879 0.00182 0.00011 0.00000 0.00054 53 20% 0.02 2.47361 0.95585 0.92538 0.16803 -18.32700 54 50% 0.05 2.43983 0.00152 -0.03538 0.00148 -6.58278 55 200% 0.2 2.43876 0.00025 0.00163 -0.00008 0.36869 56 500% 0.5 2.43875 0.00075 0.00179 -0.00006 0.37778 57 20% 0.01 2.43881 -0.00022 0.00001 0.00001 -0.00017 58 50% 0.025 2.43880 0.00096 -0.00008 -0.00002 0.00176 59 200% 0.1 2.43880 0.00122 -0.00004 -0.00001 0.00225 60 500% 0.25 2.43881 0.00003 -0.00009 -0.00002 0.00070 Elevation Pipe Changes to 848 1 0.005 0.02 2 0.005 0.1 3 0.00005 0.01 4 0.09 0.005 5 0.1 0.02 0.1 7 0.12 0.079 8 0.01 0.007 15 0.00001 0.05 974 11 0.0005 0.01 12 0.00001 0.05 13 0.00007 0.01 Transient State results NO CONVERGENCE IN CONDUIT Changes to Transient State (Refernce Scenario) % S c e n a ri o 14 0.12 0.1 9 0.05 0.01 10 0.09 0.05 890 6 0.1 111 Table 5.7. Scenarios with the highest changes of pipes water discharge relative to the reference scenario S ce n a ri o s In d ex E le v a ti o n Reference pipe Pipe changes Differences to transient state (RS) No. Hydraulic conductivit y (m/day) Diamete r (m) Diamete r % Diameter (m) Pipe water discharge % 19 848m 5 0.1 0.02 200% 0.04 23.62 20 500% 0.1 63.46 21 890m 6 0.1 0.1 20% 0.02 -19.12 25 7 0.12 0.079 20% 0.0158 -20.81 26 50% 0.0395 -12.78 28 500% 0.395 NO CONVERGENCE IN CONDUIT 36 9 0.05 0.01 500% 0.05 14.90 37 10 0.09 0.05 20% 0.01 -18.31 38 50% 0.025 -14.97 39 200% 0.1 22.48 40 500% 0.25 25.34 53 974m 14 0.12 0.1 20% 0.02 -18.33 Figure 5.10. Pipe water discharge changes for the scenarios with highest changes as listed in Table 5.7. 112 5.4.2.2 Drainage results Changes of pipe diameter D cause no tangible differences on the amount of drainage by the pipes. Except for three scenarios listed in Table 5.8, drainage changes of all other 57 scenarios are less than 0.03% relative to that of the RS. These three scenarios correspond to the highest changes scenarios of water discharge through the CFP- pipes as well. Table 5.8. Scenarios with highest drainage differences, relative to that of the RS. Scenario Difference to Reference Scenario (%) No. (from Table 5.6) Name Drainage Pipe water discharge 19 848_Pipes_5_200% 0.094 23.62 20 848_Pipes_5_500% 0.148 63.46 53 974_Pipes_14_20% 0.168 -18.33 5.4.2.3 Discharges at downstream boundaries (Constant Head OUT) The total water discharge through the constant head boundary at the downstream changes mainly only by adjustments of the diameters in pipes 9, 10 and 14 (see Table 5.6). For adjustments in pipes 9 and 10, increased downstream discharge is observed when their diameters are increased above those of the RS, whereas for pipe 14 such an increase occurs, when its pipe diameter is decreased to 20% of that of the RS. 5.4.2.4 Storage The saved and discharged amount of water in the entire system is marked in the model budget results as output storage. The highest changes of storage occur for diameter adjustments in pipes 9 and 14. For pipe 9 the storage shows big changes for all four corresponding diameter scenarios, particularly, for scenario 35, and which highlights the impact of fault F4 UP. 5.4.2.5 Discussion of pipe diameter variation scenarios Almost all scenarios which demonstrate extreme changes of pipe water discharges show also high changes of the other parameters, like drain and discharge at the downstream boundaries. There are, of course, other scenarios where one or more other parameters reach the highest value, like scenario 9, for instance, but it is not a common trend. Furthermore, the increasing of the pipe diameter has a limit with regard to the model’s numerical stability, as, for example, for scenario 28, when increasing the diameter D by 500% to 0.395 m, the model does not converge anymore. In this case, the hydraulic conductivity K of the pipe and the surrounding zones should be considered. 113 The results of Table 5.6 indicate further that for scenarios where D and K are higher changes in the pipe discharges and other budget parameters considered are more significant: At elevation 848, for instance, pipes 2 and 5 cross almost the same geological zones of the Asmari 1 formation, but due to the higher K- value of pipe 5 (0.1 m/day) than that of pipe 2 (0.005 m/day), the pipe discharge change is higher; although pipe 5’s diameter D is five times lower than that of pipe 2. The most critical pipes at elevation 890 m are pipes 7, 9 and 10 which all are located downstream of the reservoir and the grout curtain, on the left abutment of the dam body. Pipe 7, in comparison to the other pipes, has a high value of both hydraulic conductivity (0,12 m/day) and diameter (0,079 m) which manifests itself in the remarkable changes of pipe discharge in scenarios 25 and 26, when decreasing the pipe diameter to 20% and 50% of the original (RS) size. However, the most significant changes happen for pipes 9 and 10. Pipe 9 represents an upper section of the same fault (F4 Up) as pipe 5 a lower section at elevation 848 m, but with a lower hydraulic conductivity. Therefore, its effects on water discharge are lower than those of the pipe 5 at elevation 848 m, meanwhile, pipe 10 at elevation 890 m is located a little further away from pipe 9 in the same geological subunits of the Asmari 1 formation but, because of a bigger diameter and hydraulic conductivity, it leads to higher proportional changes in its water discharge, when altering its pipe diameter (see scenarios 37, 38, 39 and 40). The reason of the big amount of water budget changes when changing the diameter of the pipes making up the two F4 Up and F4 Down, is that these two long faults intersect the critical geological subunits of Asmari AS1-g and, to some extent AS1-I, both having a big karst potential (see Table 5.1). This fact has already been noted in the study phase, i.e. from the final engineering geology report of the consulting engineer Mahab Ghods (2002): “The faults F1 (not considered here, as it is locate far upstream for the dam body), F2, F3 and F4 on the left bank and F15 at the right bank of the dam should be considered in terms of water escape.” On the left bank, consequently, the most worrying fault is F4, because, in certain places, it is crossing the Asmari subunits AS1-g and As1-i, which causes water flow through the conduit paths. On the other hand, although fault F7 has contact with almost the same subunits of Amsari1 as F4, due its low wall conductivity K, its influence on pipe discharge is negligible . The drainage curtain is located at the bottom layer of the model (layer 4). Therefore, it is expected to influence significantly the pipes’ effects at elevation 848 m there. But, in fact, this happens only for pipe 5, which is again explainable with the high K- value of this pipe, in 114 comparison to pipe 2 at the same elevation and around the same position, crossing the drainage curtain on the left bank. Generally, the drainage curtain layout has no direct influence on the pipe discharge changes, neither in the corresponding layer, nor in the other top layers. Based on these above-mentioned arguments, the most determinative parameter of water flow and discharge in the modeled dual porosity formation is not the hydraulic conductivity of the porous matrix media, but the hydraulic conductivity and the size of the diameter of the pipes representing a particular fault. 5.4.3 Scenarios with hypothetical pipe locations to represent possible critical faults The previous results of the scenarios with changing pipe diameters reveal the most critical possible paths of water escape through the existing faults. To explore other probable water seepage paths, the results of the geological studies of the dam area have been perused in more detail with regard to the existence or development of other possible fault zones. These will be built into the previously calibrated CFP numerical model by means of extra hypothetical pipes to study if they lead to critical water seepage that could compromise the reliability of the dam. 5.4.3.1 Model- and pipe specifications Hence, in this section new pipes, with high values of diameter and hydraulic conductivity, are implemented into the CFP-model. The layout of these new pipes is extracted not only from the information on the remaining faults of three confirmed layers, but also from the local common borders between the various subunits of the Asmari1 formation. These remaining faults are those beyond the 15 faults already retained in the calibration of the CFP-model in Chapter 4 and were not considered there, since they had no positive impact on the calibration of the model and did not cause a decrease of the residuals and/or the rmsr between observed- and simulated heads. However, as they are indeed existing faults, they should be investigated further, not to improve the entire model, but to find out possible critical areas at the dam site. The subunits’ common borders are also considered as pipes, because movements of geological plates between these subunits could provide potential paths for water escape. As examining these hypotheses in practice in the dam area is very expensive and time-consuming, using the present numerical model will allow to test if the supposed common borders could cause water escape. Applying a try-and-error method, 10 new paths, defined by, namely, pipes 16 to 25 (following already adjusted pipes 1 to 15 in the previously analyzed transient state model), are assigned to the CFP- model as pipes in the following order, with specifications as listed in Table 5.9: 115 Table 5.9. Geological specifications of hypothetical pipes in the three layers (bottom to top) Layer/Top elevation (m) Pipe number Representative fault Common border subunits of Asmari 1 Layer 4/ 848m 16 F4 Down - Layer 3/ 890m 17 F3 - 18 - AS1-g , AS1-h 19 - AS1-h , AS1-i 20 - AS1-i , AS1-j 21 - AS1-e , AS1-f 22 - AS1-d , AS1-e 23 - AS1-f , AS1-g Layer 2/ 974m 24 F3 - 25 F4 Down - Pipe 16 in layer 4, pipe 17 in layer 3 and pipes 24 and 25 in layer 2 are representatives of faults on the left abutment. Pipes 18 to 23 are located at the common borders of the different Asmari 1 subunits. Except for pipes 24 and 25 on the right abutment, all other pipes are set at the left abutment (see Figure 5.11). For almost all pipes three scenarios are defined: 1) the diameter of the pipe is assigned as D=0.1 m and the hydraulic conductivity K to 0.1 m/day, 2) the diameter is doubled to D=0.2 m, using the previous K, and 3) the primary value of D is kept again (=0.1m), but K is increased by 100% to 0.2 m/day. The above specifications are based on pipe reference values found by trial and error in the adjusted pipe parameter simulations of the previous section. As pipe 17 is the continuation of an existing pipe 7 of the reference scenario, the primary values of that pipe are the same as that of pipe 7 (D=0.079 m, K=0.12m/day). 5.4.3.2 Results and discussions Table 5.10 summarizes the model results for all above-mentioned scenarios, including the changes of cumulative water discharges (at the end of the 1727-day simulation period) at the downstream boundary in comparison to the reference scenario (RS). The hypothetical pipes with the highest changes of water discharges are illustrated in the bar diagram of Figure 5.12. The results summarized in Table 5.10 show the importance of the primary choice of the assigned pipe parameters, namely, diameter D and hydraulic conductivity K on the total water discharge across the downstream boundary of the model. Initially, D= 0.1 m and K = 0.1 m/day are obtained for pipes 1 to 15 from the calibration of the transient CPF- model and the results of the pipe diameter variations in Section 5.4.1. In all, but some scenarios of pipes 18 and 24, there are very noticeable effects of the adjusted hypothetical pipes on the rmsr and the outflow water discharge. 116 Figure 5.11 Layout of hypothetical- and adjusted pipes in a) layer 2, b) layer 3 and c) layer 4; and d) layout of all hypothetical pipes in plan view. The barplots of the water discharge changes of the most influential hypothetical pipes (Figure 5.12), with their locations depicted in Figure 5.11, indicate that the left abutment is the critical side of the dam area, as the water escape changes through the downstream boundaries are the largest when including one of these pipes into the CFP- model. 117 Table 5.10. Changes of water escape in comparison to RS for the three hypothetical pipes scenarios Layer/ Top Elevation Pipe number Scenario number Diameter Hydraulic Conductivity Downstream water discharge Changes rmsr m m/day m3 m3 % m Transient State (Reference Scenario) results 90,914,165 - - 2.439 Layer 4/ 848m 16 61 0.1 0.1 90,906,834 -7,331 -0.008 2.444 62 0.2 0.1 90,926,736 12,571 0.014 2.450 63 0.1 0.2 90,844,538 -69,627 -0.077 2.523 Layer 3/ 890m 17 64 0.079 0.12 90,954,163 39,998 0.044 2.441 65 0.158 0.12 NO CONVERGENCE - 66 0.079 0.24 90,985,399 71,234 0.078 2.442 18 67 0.1 0.1 90,918,618 4,453 0.005 2.436 68 0.2 0.1 90,891,269 -22,896 -0.025 2.470 69 0.1 0.2 91,747,556 833,391 0.917 2.433 19 70 0.1 0.1 91,749,037 834,872 0.918 2.471 71 0.2 0.1 91,749,441 835,276 0.919 2.471 72 0.1 0.2 91,764,067 849,902 0.935 2.469 20 73 0.1 0.1 91,760,820 846,655 0.931 2.473 74 0.2 0.1 91,761,434 847,269 0.932 2.473 75 0.1 0.2 91,789,349 875,184 0.963 2.472 21 76 0.1 0.1 91,405,853 491,688 0.541 2.424 77 0.2 0.1 91,407,534 493,369 0.543 2.424 78 0.1 0.2 91,435,256 521,091 0.573 2.472 22 79 0.1 0.1 91,429,223 515,058 0.567 2.444 80 0.2 0.1 91,430,210 516,045 0.568 2.444 81 0.1 0.2 91,466,276 552,111 0.607 2.444 23 82 0.1 0.1 91,368,600 454,435 0.500 2.466 83 0.2 0.1 91,368,560 454,395 0.500 2.466 84 0.1 0.2 91,366,829 452,664 0.498 2.486 Layer 2/ 974m 24 85 0.1 0.1 90,896,235 -17,930 0.020 2.439 25 86 0.1 0.1 91,365,734 451,569 0.497 2.445 87 0.2 0.1 91,365,723 451,558 0.497 2.445 88 0.1 0.2 91,365,637 451,472 0.497 2.446 Layer 3, ranging from elevation 848 to 890 m.a.s.l, is the most critical layer of the model, as it intersected by most of the active hypothetical pipes. The top layer 1, between 974 and 1032 m.a.s.l is unconfined and most cells downstream of the reservoir are dry. Although the second top layer 2 (890 to 974 m.a.s.l) is confined, in many of its cells downstream of the reservoir, the water levels are close to the elevation of the layer’s bottom. 118 Figure 5.12 Water leakage changes for hypothetical pipes with the largest effects (see Table 5.10) Three critical pipes, namely, pipes 18, 19 and 20 in layer lead to high changes of water escape of more than 800,000 cubic meters at the end of the simulation period which, considering that the total amount of the RS’ water discharge is about 91 million cubic meters, means that these pipes lead to an increase of discharge through the downstream boundary of nearly 1% of the total discharge of the RS. In fact, these three pipes are all located at the common border between the different Asmari 1 subunits. The most critical scenario of each of these three pipes is the third one, where the hydraulic conductivity K has been doubled, and which so explains the importance of the latter in defining the flow exchange between the matrix medium and the discrete conduit. But it does not necessarily mean that increasing K will lead to a proportional rise of the water escape. Moreover, pipe 18 shows no noticeable water escape for its first- (primary) and second (doubling of the pipe diameter D) scenarios, whereas, for pipes 19 and 20, high changes of water discharge are obtained for all three scenarios, with the third one (doubling of K) leading to the highest discharge changes of all hypothetical pipe scenarios. Similarly, for the pipe 24- scenarios the model fit measured by the rmsr indicated a large amount of insensitivity to the inclusion of this pipe, so no other than the first baseline scenario were executed. The reason for the high amounts of water escape when introducing these three pipes in layer 3, which have a common border with subunits of AS1-g for pipe 18 and AS1-i for pipes 19 and 20, is, in addition to the wall-conductivity of the pipes, also the high Darcy hydraulic conductivity of the surrounding matrix medium. These two subunits have almost the same geological specifications, particularly, a high karstification potential and so a very high 119 permeability (Mahab Ghods report, 2002). Furthermore, they have an intersection with the most critical pipe 9 (fault F4-Up) in layer 3 of the calibrated transient state CFP-model (see Figure 5.11), which will further promote water flow through these paths. None of the scenarios of adjusted hypothetical pipes at the locations of the exiting faults (pipes 16, 17 and 24), result in a high amount of water escape (see Table 5.10). Therefore, the faults simulated by these pipes should not be of concern in terms of significant water escape (seepage) in the Karun 4 dam area. This means further that the grout curtain at the positions of these pipes/faults works satisfactorily as a sealing obstacle to water seepage from the dam/reservoir. All these results taken together confirm the impacts of the geology, as analyzed and reported in the previous study of the authors (Hosseiny Sohi et al., 2017a) on the outcome of the MODFLOW-CPF model. 5.4.4 Scenarios with connected pipes (continuous faults) through the grout curtain 5.4.4.1 Rational and approach To evaluate the impact of the earlier recognized faults with their pipe-representations (see Table 5.3 and Figure 5.7, Section 5.3.2.2) on the grout curtain and its water sealing capacity, some scenarios are defined by joining those pair-pipes at the up- and downstream of the reservoir which pass through the grout curtain. This corresponds to assuming that the sealing of the latter is not good enough to withhold water in these parts so that the opened conduits will cause downstream water escape. The realization of these scenarios is effectuated by connecting pipes 1 and 5 in layer 4, pipes 8 and 9 in layer 3, and pipes 12 and 11 in layer 2, upstream to downstream of the reservoir across the grout curtain. All these pipes/faults are located in the left abutment (see Figure 5.13). To that regard, the connection pipes 1-5 and pipes 8-9 are representatives of the continuous fault F4-Up in layer 4 and layer 3, respectively, whereas pipes 12-11 are simulating fault F2 in layer 2. In the final calibrated transient state (reference scenario, RS) model, these pipes have different diameters and hydraulic conductivities (see Table 5.3). To establish the similar characteristics, two scenarios are set up to each of these connected pipes: assuming, firstly, the pipe diameter D and hydraulic conductivity K of the upstream pipe; and, secondly, those of the downstream pipe, as representative for the whole connected pipe. These variants result in a total of six scenario model runs. 120 Figure 5.13 Layout of three connected pipes in three layers in the left abutment In addition, two universal scenarios with all three connected pipes together in the model are executed: one with the minimum- and another one with the maximum of the three pipe parameter values. And finally, the last scenario is designed based on the results of Section 5.4.3, where it has been indicated that the highest impact of the hypothetical pipe on the model’s downstream water discharge occurs for a pipe diameter D=0.1 m and a hydraulic conductivity K =0.2 m/day. 5.4.4.2 Results For all of the 9 scenarios mentioned, the water discharge at the downstream boundaries of the model is evaluated again and compared with that of the RS. The results of these scenarios (#89-97) are presented in Table 5.11. The table shows that, compared with the volume of the water escape of the RS, none of the connected pipes in the individual scenarios 89-94 leads to a noticeable increase of water escape, as these are all lower than 150,000 m3, while the corresponding values for the scenarios of the combined connection of the three pipes across the grout curtain, i.e. pipes 1-5, 8-9 and 12-11 in scenarios 95 and 96 amounts to more than 800’000 m3. Nevertheless, in both scenarios, the water discharge amounts are almost the same as those obtained for scenarios 69-75 of the hypothetical pipes 19, 18 and 20 in Table 5.10. 121 Table 5.11. Results of different scenarios of the connected pipes at the up- and downstream sections of the reservoir across the grout curtain. S ce n a ri o Description Pipe Water discharge (m3) Number Hydraulic conductivity (m/day) Diameter (m) RS Difference relative to RS 90,914,165 89 Connected_Pipes_1_848_5 1 0.005 0.02 90,927,225 13,060 90 Connected_Pipes_5_848_1 5 0.1 0.02 90,933,414 19,249 91 Connected_Pipes_8_890_9 8 0.01 0.007 91,063,511 149,346 92 Connected_Pipes_9_890_8 9 0.05 0.01 90,920,084 5,919 93 Connected_Pipes_11_974_12 11 0.0005 0.01 90,914,793 628 94 Connected_Pipes_12_974_11 12 0.00001 0.05 90,914,071 -94 95 Three connected pipes min Connected_Pipes_1_848_5 0.005 0.02 91,746,353 832,188 Connected_Pipes_8_890_9 0.01 0.007 Connected_Pipes_11_974_12 0.0005 0.01 96 Three connected pipes max Connected_Pipes_5_848_1 0.1 0.02 91,755,602 841,437 Connected_Pipes_9_890_8 0.05 0.01 Connected_Pipes_12_974_11 0.00001 0.05 97 Three connected pipes hypothetical Connected_Pipes_5_848_1 0.2 0.1 93,989,871 3,075,706 Connected_Pipes_9_890_8 0.2 0.1 Connected_Pipes_12_974_11 0.2 0.1 The above results explain that, although the connected pipes represent the most concerning fault F4-Up (see Figure 4.7) in the four scenarios (89 to 92), they have individually not a significant impact on the amount of water escape out of the groundwater flow system. Hence it can be deduced that the installed grout curtain underneath the dam works efficiently and the risk of undesirable water leakage is not high. On the other hand, for the more critical pipe-connection cases, where all three pipes at the up- and downstream of the grout curtain act in parallel (scenarios 95 and 96), the additional water discharge volume across the downstream boundary of the model will be significant, albeit still less than one percent of the entire of the total water discharge in the model. Finally, regarding the last, extreme scenario 97, which truly leads to a remarkable increase of the amount of water leakage out of the system, this scenario may not be considered as very realistic for representing faults cutting through the grout curtain, as the probability of having such high values for both the hydraulic conductivity and the pipe diameter of all three pipes (K=0.2 m/day and D=0.1 m), appears to be rather low. 122 5.5. Conclusions Different scenarios of possible water leakage through the karstic dam foundation and abutments of Karun 4 dam have been investigated in this study. To that avail, the numerical groundwater flow model MODFLOW-CFP is employed that simulates laminar groundwater flow in the porous matrix medium and laminar or turbulent flow in discrete conduits representing faults/fractures/fissures which, in the present application, are mostly due to karstification processes. 15 faults around the dam site, recognized from original geological site reports, are incorporated into the model in the form of pipes of given diameter and exchange/wall hydraulic conductivity which are further fine –tuned during the trial and error calibration process. To that avail, observed data of piezometers, observation wells and drainage from the dam site is employed. This calibrated and validated transient model defines then the reference scenario (RS) that represents the current situation of the groundwater flow in the area. For a deeper understanding of the effects of the existing and other potential karst conduits on the water leakage from the dam, in the subsequent exhaustive analysis, a total of 97 conduit scenarios, divided into three categories, are simulated and the computed amount of downstream water escape compared with that of the transient reference scenario above: (1) adjusting the diameters of each of the 15 pipe representing the most relevant faults fourfold of the corresponding reference diameter of the pipe in the model result in scenarios 1-60, (2) embedding of 10 new hypothetical pipes at locations of extra faults not considered previously with three combinations of diameter and hydraulic conductivity result in another 28 scenarios (2 not considered) (scenarios 61-88), (3) connecting several adjusted pair-pipes to each other to represent continuous existing and not-yet identified faults on the up- and downstream section of the grout curtain (scenarios 89- 97). From the results of the first group of scenarios, namely, changing the diameters of the pipes, the most critical fault, the one leading to the highest amount of downstream water escape out of the system could be identified, which turned out to be fault F4-Up. The reason for this potentially detrimental behavior of this fault is that it intersects both AS1-g and AS1-i subunits of the Asmari 1 formation with its large degree of karstification (see Table 5.1), confirming results of the previous geological investigations (Hosseiny Sohi et al., 2017a). 123 The results of the second group of scenarios, with new hypothetical pipes introduced at locations of other potential faults, demonstrate water discharge increases of about 1% relative to the reference scenario for faults at the border of the adjacent subunits As1-g and As1-i of the Asmari-1 formation. The results of the connected-pipes scenarios (group three) serve mainly to understand the effectiveness of the grout curtain in blocking downstream water seepage, in particular, as there is the earlier recognized critical fault F4-Up cutting through it. Thus, the corresponding scenario results indicate that, while connecting pipes individually across the grout curtain to mimic potential intersecting faults, has practically no noticeable, additional effect on the water discharge out of the groundwater flow system, increases of less than 1% are obtained only when assuming all three pipes acting in parallel, which appears to be unrealistic. Hence it can be deduced that the installed grout curtain underneath the dam works efficiently, and the risk of undesirable water leakage is relatively low. This is the more so, since the reservoir reached its maximum water level in several stress periods during the modelling period of nearly 5 years, which means that the above scenario simulations do in fact include the most extreme boundary conditions for driving water seepage in the model domain. In summary, the results of the application of the novel MODFLOW-CFP model to the Karun 4 dam site lead to the following final conclusions: - Possible water leakage paths in the dam area are identified and their impacts on the amount of downstream seepage investigated. - The left abutment of the dam site has a higher potential of water leakage, due to the presence of some pronounced karstic fault. - The aquifer layer between 848 to 890 m.a.s.l (layer 3 of the model) is the most critical geologic unit, not only in the left abutment, but also in whole model area due to its high degree of karstification. - Under the present reservoir conditions, the grout curtain appears to serve its purpose to impede water leakage, in spite of the presence of some critical faults. In spite of the above, seemingly benevolent water leakage situation underneath the Karun 4 dam at present, further ongoing monitoring of the dam area is strongly recommended, namely, of the various critical karstic faults, whose effects on downstream water seepage have been clearly demonstrated here employing the novel MODFLOW- CFP- dual porosity model. 124 References Anderson MP, Woessner WW & Hunt RJ (2015). Applied groundwater modeling: Simulation of flow and advective transport, 2nd Edition. Assari A & Mohammadi Z (2017). Assessing flow paths in a karst aquifer based on multiple dye tracing tests using stochastic simulation and the MODFLOW-CFP code, Hydrogeology Journal, Vol 25, No 6, pp 1679-1702. Ashok KC & Sophocleous MA (2009). Recent MODFLOW developments for groundwater modeling: Kansas Geological Survey Open-file Report 2009-4, 139 p. Atkinson TC (1977). Diffuse flow and conduit flow in a limestone terrain in the Mendip Hills, Somerset, J. Hydro., 19, 323– 349. Barenblatt GK, Zheltov IP & Kochina N (1960). Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks, Prikl. Mat. Mekh., 24, 5, 852– 864. Brown GO (2003). The History of the Darcy-Weisbach Equation for Pipe Flow Resistance, In Rogers, J. R.; Fredrich, A. J. Environmental and Water Resources History. American Society of Civil Engineers. pp. 34–43. Clemens T, Huckinghaus D, Sauter M, Liedl R & Teusch G (1996). A combined continuum and discrete network reactive transport model for simulation of karst development. Proceedings of the ModelCARE 96 Conference held at Golden, Colorado, September 1996. IAHS Publ. no. 237, 306- 318. Cornaton F & Perrochet P (2002). Analytical 1-D dual-porosity equivalent solutions to 3-D discrete single continuum models. Application to karstic spring hydrograph modeling, J. Hydrology, 262, 165-176. Darcy H (1856). Les Fontances publiques de la ville de Dijon: Paris, Victor Dalmont. Ford D & Williams P (2007). Karst Hydrogeology and Geomorphology. Wiley Ltd, England. Gallegos JJ, Hu BX & Davis H (2013). Simulating flow in karst aquifers at laboratory and sub-regional scales using MODFLOW-CFP: Hydrogeology Journal, v. 21, p. 1749-1760. Ghasemizadeh R, Hellweger F, Butscher CH, Padilla I, Vesper D, Field M & Alshawabkeh A (2012). Review: Groundwater flow and transport modeling of karst aquifers, with particular reference to the North Coast Limestone aquifer system of Puerto Rico. Gunn J (1985). a conceptual model for conduit flow dominated karst aquifers, Karst Water Resources (Proceedings of the Ankara - Antalya Symposium, July 1985). IAHSPubl. no. 161. Haghnejad A, Ahangari K & Noorzad A (2014). Investigation on various relations between uniaxial compressive strength, elasticity and deformation modulus of Asmari formation in Iran. Arab. Journ. Science and Engineer., 39, 4, 2677-2682. Halihan T & Wicks CM (1998). Modeling of storm responses in conduit flow aquifers with reservoirs, J. Hydrol., 208, 82-91. 125 Harbaugh AW (2005) MODFLOW-2005, the U.S. Geological Survey modular groundwater model— The ground- water flow process: U.S. Geological Survey Techniques and Methods 6-A16. Hartmann A, Goldscheider N, Wagener T, Lange J & Weiler M (2014). Karst water resources in a changing world: review of hydrological modeling approaches. Rev Geophys 52:218–242. Hill ME, Martin A & Stewart MT (2008). Performance Evaluation of the Modflow-2005 Conduit Flow Process Applied to a Karst Aquifer Underlying West-Central Florida, in E. L. Kuniansky, editor, U.S. Geological Survey Karst Interest Group Proceedings, Bowling Green, Kentucky, May 27-29, 2008: U.S. Geological Survey Scientific Investigations Report 2008-5023, p.93-98. Hosseiny Sohi SM, Koch M & Ashjari J (2014). Construction Effects of the Karun 4 Dam, Iran, on the Groundwater in the Adjacent Karstic Aquifer. ICHE 2014, Hamburg - Lehfeldt & Kopmann (eds). Hosseiny Sohi SM, Koch M & Ashjari J (2017a). Evaluating permeability and groutability at Karun 4 dam Iran using Lugeon values and grout Take. Symposium Proceeding of 85th Annual Meeting of International Commission on Large Dams, Prague, July 2017. Hosseiny Sohi SM, Koch M & Ashjari J (2017b). Monitoring of the seepage controlling system of the Karun 4 dam in Iran. Symposium Proceeding of 85th Annual Meeting of International Commission on Large Dams, Prague, July 2017. Houlsby AC (1990). Construction and design of cement grouting. A guide to grouting in rock foundations. New York, Wiley (Wiley series of practical construction guides). Huntoon P W (1995). Is it appropriate to apply porous media groundwater circulation models to karstic aquifers El-Kadi, A. I., ed., Groundwater models for resources analysis and management, 1994 Pacific Northwest/Oceania Conference, Honolulu, HI, p. 339-358. Jaeger Ch (1979). Rock mechanics and engineering. Cambridge University Press. Jeannin P-Y (2001). Modeling flow in phreatic and epiphreatic karst conduits in the Holloch Cave (Muotatal, Switzerland). Water Resour. Res., 37, 191-200. Keeler RR & Zhang YK (1997). Modeling of groundwater flow in a fractured-karst aquifer in the Big Springs Basin, Iowa. GSA Abs. with Programs, 29, 4, 25. Kiraly l (1998). Modeling karst aquifers by the combined discrete channel and continuum approach. Bulletin du Centre d’hydrogeologie, Neuchatel, 16, 77-98. Knill, JL (1972). Assessment of reservoir feasibility. Quart. J. Eng. Geol, 4, 355-372. Kuniansky EL & Shoemaker WB (2009). Conduit Flow Process (CFP) for MODFLOW- 2005: in Proceedings 15th International Congress of Speleology, Kerrvile, Texas, USA, July 19- 26, 2009; William B. White, ed.: Greyhound Press, Huntsville, Alabama: pp. 1568-1574. Mahab Ghods Consulting Engineering Co. (2002). Final engineering geological report (phase 2), Karun 4 dam project. Mahab Ghods Consulting Engineering Co. (2010). Final excavation and injection report of the Karun 4 dam grout curtain. Mahab Ghods Consulting Engineering Co. (2015). Monitoring report of Karun 4 dam. 126 Milanović PT (2004). Water resources engineering in karst. CRC Press, Boca Raton. Mohrlok U & Teutsch G (1997). Double continuum porous equivalent (DCPE) versus discrete modeling in karst terranes. In: Gunay G, Johnson AI (eds) Karst waters and environmental impacts. In: “Proceedings 5th Int Symp Field Sem Karst Waters and Environmental Impacts, Turkey. Balkema, Rotterdam, pp. 319–326. Neuman SP (2005). Trends, prospects and challenges in quantifying flow and transport through fractured rocks. Hydrogeology Journal 13 (1), 124 - 147. Quinlan JF & Ewers RO (1985). Groundwater flow in limestone terrains: Strategy rationale and procedure for reliable and efficient monitoring of groundwater quality in karst areas. In: Proceed. 5th national symposium and exposition on Aquifer Restoration and Groundwater Monitoring, The Fawcett Center, Columbus, Ohio. pp 197-234. Reimann T, Geyer T, Shoemaker WB, Liedl R & Sauter M (2011a). Effects of dynamically variable saturation and matrix-conduit coupling of flow in karst aquifers, Water Resources Research, v. 47, W11503, doi:10.1029/2011WR010446. Reimann T, Birk S, Rehr C & Shoemaker W.B (2011b). Modifications to the Conduit Flow Process Mode 2 for MODFLOW-2005: Ground Water, doi: 10.1111/j.1745-6584.2011.00805.x Reimann T & Hill ME (2009). MODFLOW-CFP: A New Conduit Flow Process for MODFLOW-2005: Ground Water, v. 47, p. 321-325. Sadd Tunnel Pars Company (2015). Evaluation of Grouting Curtain of Karun 4 dam. Saller SP, Ronayne MJ & Long AJ (2013). Comparison of a karst groundwater model with and without discrete conduit flow: Hydrogeology Journal, v. 21, no. , p. 1555-1566. Sauter M (1992). Quantification and forecasting of regional groundwater flow and transport in a karst aquifer (Gallusquelle,Malm,SW Germany): Tübinger Geowissenschaftliche Arbeiten, Ser.C, 13, 150p. Shoemaker W B, Kuniansky EL, Birk S, Bauer S & Swain ED (2008). Documentation of a conduit flow process (CFP) for MODFLOW-2005, U.S. Geol. Surv. Tech. Methods, Book 6, Chap. A24. Uromeihy A (2000). The Lar Dam; an example of infrastructural development in a geologically active karstic region. Journal of Asian Earth Sciences, Volume 18, Issue 1, pp 25–31. Weisbach J (1845). Lehrbuch der Ingenieur und Maschinen Mechanik, Vol. 1. Theoretische Mechanik, Vieweg und Sohn, Braunschweig. 535 pages (in German). White WB (1999). Conceptual models for karstic aquifers. In: Karst Modeling, Karst Waters Inst. Spec. Publ., Vol. 5, A. N. Palmer, M. V. Palmer and I. D. Sasowsky (eds), pp. 1-16, Karst Waters Inst., Charles Town, W. Va. Winston RB (2009). ModelMuse: a graphical user interface for MODFLOW-2005 and PHAS T. US Geological Survey, Reston, VA. Worthington SRH, Ford DC & Beddows PA (2000). Porosity and permeability enhancement in unconfined carbonate aquifers as a result of solution, In “Speleogenesis: Evolution of Karst Aquifers” A. B. Klimchouk et al. (eds)., pp. 463-472, Natl. Speleology Soc., Huntsville, Al. 127 Yurtsever Y & Payne R.B (1986). Mathematical models based on compartmental simulation approach for quantitative interpretation of tracer data in hydrological systems, Proc. 5th Intl. Symposium of Underground Water Tracing, IGME, Athens, Greece, pp. 341-353. Zhang Y-K, Bai E-W, Libra R, Rowden R & Liu H (1996). Simulation of spring discharge from a limestone aquifer in Iowa, Hydrogeol. J.., 4, 41-54.