Seepage analysis of Upper Gotvand Dam concerning gypsum karstification (2D and 3D approaches)

Abstract

Upper Gotvand Dam is constructed on the Karun River at the south west of Iran. In this paper, 2D and 3D models of the dam together with the foundation and abutments were established, and several seepage analyses were carried out. Then, the gypsum veins that are scattered throughout the foundation ground were included in the models, and the seepage pattern, considering the dissolution law of gypsum, was analyzed. It was disclosed that the discharge fluxes obtained from 2D and 3D analyses are not similar, and the discharge flux in 3D model is about four times that of the 2D model. Also, the 3D model locates the phreatic surface somewhat higher than the 2D model. This means that the 2D model estimates lower pore water pressure pattern in comparison with the 3D model. These may be attributed to the fact that with 2D model the lateral components of vectors of seepage velocity are ignored. In the current case, the rate of increase of discharge flux due to dissolution of gypsum veins was obtained to be a third-order function of the aperture width. In spite of the fact that the grout curtain is designed to be about 170 m deep, however, complete dissolve of gypsum will severely increase the discharge flux through the foundation ground.

Keywords

Upper Gotvand Dam / seepage analysis / gypsum veins / 3D model / discharge flux

Introduction

It is estimated that gypsum or anhydrite deposits underlie approximately 25% of the land surface. Only 10 % of these deposits outcrop. At these outcrops, or where gypsum or anhydrite strata occur in depths of a few hundred meters, gypsum karst has evolved. Therefore, extensive areas of gypsum karst exist worldwide [1,2]. Some of the problems caused by the dissolution of gypsum and anhydrite in a dam construction site would be such as gradual increase of seepage rate through dam abutments and foundation ground, dam breakage because of intensive leakage through foundations, making the dam reservoir unusable because of intensive water escape through large leakage paths, strength reduction in foundation ground because of the gypsum dissolution, and concrete structures destruction because of the sulfated water caused by dissolution of gypsum and anhydrite [3]. Saint Francis and San Fernando in California State of the USA are two examples of the damaged dams due to the gypsum and anhydrite wash away. One of the solutions to prevent such problems is to deploy seepage analysis on a realistic model of embankment together with the foundation and abutments. Although the three-dimensional analysis is an appropriate tool for these cases, however, its application to earth dams is rather complicated. According to engineering practice, in order to estimate the seepage rate usually, one or more critical sections are employed and analyzed two dimensionally. Generally, the two dimensional analysis, especially in narrow valleys or valleys with varying profile, is erratic. The main reason is that the lateral component of seepage velocity that is ignored in 2D approach may be quite considerable. Although the Massingir Dam has a uniform section in 3 km length, with a simple three-dimensional seepage analysis, more realistic results were obtained compared with two dimensional analysis, and the importance of three-dimensional seepage analysis was indicated [4].

Many studies about seepage analysis through dam abutment or foundation are performed. Uromeihy et al. [5] evaluated the potential of water seepage at the Chapar-Abad Dam by the study of joint systems of the rock units, the use of numerical analysis to simulate water flow in the ground, and by conducting in-situ tests to estimate the permeability’s values. Ghobadi et al. [6] investigated the seepage phenomena through the right abutment of Shahid Abbaspour Dam. Ahmed [7] investigated the unconfined flow through dams and modeled hydraulic conductivity as spatially random field following lognormal distribution.

In this paper, 2D and 3D models of the dam together with the foundation and abutments were established, and several seepage analyses were carried out. Then, the gypsum veins that are scattered throughout the foundation ground were included in the models, and the seepage pattern, considering the dissolution law of gypsum, was analyzed.

Three-dimensional flow

The three-dimensional flow equation is derived, which is called coupled equation of flow. For an unsaturated soil having heterogeneous anisotropic conditions, the coefficient of permeability at a point varies in x, y, and z directions. However, the permeability variations in the three dimensions are assumed to be governed by the same permeability function. The continuity for three-dimensional steady-state flow can be satisfied as follows [8]:

(1)\left({\displaystyle \frac{\partial v_{\mathrm{w}x}}{\partial x}}+{\displaystyle \frac{\partial v_{\mathrm{w}y}}{\partial y}}+{\displaystyle \frac{\partial v_{\mathrm{w}z}}{\partial z}}\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z=\mathrm{0},

where ν_wx, ν_wy, and ν_wz are water flow rate across a unit area of the soil in x, y, and z directions, respectively. Referring to Eq. (1), the governing differential equation for the steady-state seepage analysis in 3D space may be derived as

(2){\displaystyle \frac{\partial }{\partial x}}\left(k_x{\displaystyle \frac{\partial H}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(k_y{\displaystyle \frac{\partial H}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left(k_z{\displaystyle \frac{\partial H}{\partial z}}\right)+Q=\mathrm{0},

where k_x, k_y, and k_z are the coefficients of permeability in x, y, and z directions, respectively; H is the total head; and Q is the flux at the model boundaries. Using the Galerkin principle of weighted residuals, the finite element formulation for steady-state seepage in three dimensions is derived as [9]

(3){\displaystyle {\int }_v{B}^{\mathrm{T}}CB\mathrm{d}v}H={\displaystyle \int Aq{N}^{\mathrm{T}}N\mathrm{d}A},

where BoldItalic and BoldItalic are gradient matrix and element hydraulic permeability matrix, respectively. BoldItalic is the vector of nodal heads, and A is the area of the face of the element. q is the unit flux across the faces of an element, and BoldItalic shows the vector of interpolating functions.

Permeability function

There are several geotechnical and geo-environmental problems involve water flow through soils. However, the coefficient of permeability is the most important parameter that dominates the water flow pattern through soils [10]. It has been shown that in the steady-state seepage problems the results obtained using a typical permeability function may be quite close to those of exact solution [11].

Figure 1 represents the coefficient of permeability as a function of pore water pressure. Assume that for the first iteration, all elements are assigned a saturated permeability (k_a) corresponding to zero pressure.

This permeability will allow for more flow than is required and will result in a highly negative pore water pressure (point P_f). For the next iteration, the permeability will be k_f. This value does not allow for enough flow, and the computed pressures will be positive. Once again, the permeability will be set to a value that is too high, resulting in a solution that oscillates between the extremities permitted by the function [9]. However, this figure clearly indicates the importance of application of a pressure dependent permeability function for both saturated and unsaturated zones of an earth dam. In this way, such a zoning procedure will result in more accurate seepage flux than the case in which if a constant hydraulic permeability function would be used [11]. Using a constant permeability value may result in erratic results. Phreatic surface position may be unrealistic, and flow rate in unsaturated zone may be extremely high. When using a constant permeability value, water can flow through unsaturated zone like that of saturated zone [9]. Hydraulic permeability values of the embankment can be measured from tests, and there are several methods introduced in different references.

Upper Gotvand Dam project

The dam is constructed across the Karun River, in the Khuzestan Province of Iran, north of Shushtar city. It is a 178 m high earth dam with central clay core. The dam is underlain by the Bakhtiari formation (BK) overlying the Aghajari formation (AJ). While the BK formation remains in its horizontal position, the AJ formation is folded and faulted, and the bedding planes and joints are inclined with varying dip angles.

Apertures within the BK formation are usually vertical and have relatively large openings and continuity of the apertures often reaches to several meters [12]. In Fig. 2, the geological layers of Upper Gotvand Dam are depicted. The AJ-rocks contain veins of gypsum usually associated with clay stone beds. The maximum thickness of the veins is reported to be 2 cm and even wider that appear as thin films on the beddings and along joint planes.

The latter ones were formed subsequently and indicate the dissolution of gypsum and re-sedimentation in the joints. In general, gypsum is encountered below 25 m deep. Inspecting foundation borehole logs of the AJ formation at the dam site, four different depths were observed in which gypsum veins exist. These are summarized in Table1.

Modeling

The Seep3D software was employed as an effective tool for seepage analysis. Seep3D is a new software product for modeling three-dimensional seepage problems. Seep3D is formulated for conditions of constant total stress, that is, there is no loading or unloading of the soil mass.

It is assumed that the pore-air pressure remains constant at atmospheric pressure during transient processes [9]. This software uses finite element method for seepage analysis, and it is able to consider unsaturated condition. This software contains three elements, namely, hexahedron, prism, and tetrahedron, for modeling, as shown in Fig. 3(a).

In order to establish an appropriate model for seepage analysis, following steps were taken.

Geometry

In order to create the geometry of the model, topographic maps of the dam site and also section drawings at the distances of 50 m were used; also, the hexagonal element in Seep3D software was employed (see Fig. 3).

There are two boundary conditions in steady-state seepage analysis. The maximum water level in the reservoir was considered as the upstream boundary surface. For downstream, potential seepage condition was applied. In other words, all surfaces that water may seepage through were considered as potential seepage surfaces.

It is common in finite element method to model the dam with the adjacent zones, including reservoir, abutments, and foundation. Since the required extension of the boundary in the foundation ground is not known, it is necessary to carry out some boundary sensitivity analyses. However, for 3D analysis, an overall mass of 200 m in depth and 300 m in width (abutments) was obtained to be quite satisfactory.

Material properties

The embankment material and foundation ground engineering properties were introduced as follows: for the embankment materials the laboratory permeability test results at zero pressure were used to establish the general permeability functions. These functions are shown in Fig. 4.

Compacted soil behaves transversally isotropic in which, k_x = k_y ≠ k_z. Generally, with earth dams, it may be accepted that [9,13] k_x / k_y = 1 and k_x / k_z = 10.

For foundation ground, some considerable numbers of Lugeon permeability test results are available. These results were averaged and used to categorize the foundation ground layers in terms of permeability. However, it should be noted that the foundation ground was considered hydraulically isotropic, i.e., k_x = k_y = k_z.

In Table 2, results of Lugeon test results at AJ formation in foundation ground are presented. These results are from exploratory boreholes in dam foundation ground. Regarding this table, the weighed average Lugeon value is estimated to be 6 Lu, which is equal to about 9 × 10^-7 m / s.

Trial grouting in AJ formation showed that cement grout could not satisfactorily penetrate into the gypsum veins [12].

Sensitivity analyses of 2D and 3D models for the number of elements

In order to select the appropriate finite element mesh, sensitivity analyses for the number of elements were carried out for both 2D and 3D models. The results are shown in Figs. 5 and 6.

It can be seen that increasing the number of elements in the 2D model decreases the discharge flux; hence, the appropriate number of elements for 2D analysis is around 17000 elements. For 3D model, results converge at about 25000 elements, however. Therefore, in seepage analysis with 3D model, convincing results can be achieved by breaking the mesh up to 32000 elements and finer.

Modeling of gypsum veins in foundation ground

Gypsum is hardly soluble in water, but with water temperature and pressure rise, its solubility increases [14]. Water flowing through narrow fissures and fractures in soluble rock, such as limestone and gypsum, widens these by chemical dissolution. When close to saturation characteristic, nonlinear dissolution kinetics cause a feedback mechanism, giving rise to a breakthrough behavior of the flow rates.

If constant head h drives flow of aggressive water through a single fracture of aperture width a₀ and length L, initially, flow rates increase slowly, but then, suddenly at breakthrough time, T, they increase dramatically. Recent modeling approaches on two-dimensional domains of dam sites have shown that under unfavorable conditions, leakage below dam sites can increase to an unbearable extent within the lifetime of the structure [15].

In case of flowing water, dissolution phenomenon continues until the gypsum has been washed away completely, leaving wider apertures and open fractures. In order to model mode of formation and performance of these apertures, the following hypothesis was used. According to Fig. 7 and Eq. (4), equivalent permeability for the supposed rock mass can be computed as shown below [16]:

(4)k={\displaystyle \frac{\rho g}{\mathrm{12}\eta }}{\displaystyle \frac{a_{\mathrm{0}}^{\mathrm{3}}}{s}},

where ρ is density of water; η is viscosity of water and equals 0.0065 at 25°C; g is acceleration due to gravity; a₀ is aperture width; and s is fractures average spacing.

In Table 1, gypsum veins are classified into four general depths. Generally, water flow through fractures is either laminar or turbulent flow. In the early stage, there is laminar flow through apertures, and dissolution is faster at the aperture entrance. As the calcium concentration in pore water increases, the rate of dissolution of gypsum reduces, and a funnel-like conduit evolves there. This opening at the exit enhances the flow rate through the fracture, and therefore, the funnel-like opening at the entrance propagates further downstream, and also, the dissolution rates at the exit increase further. By this time, the water flow becomes turbulent, and now, the flow rate is so high that the concentration of calcium becomes close to zero, and dissolution rate tends to be even along the fracture. Therefore, as time proceeds, the funnel-like shape becomes smoothened out [16]. Figure 8 shows the variations of the AJ formation mass permeability against aperture width. This figure reveals that dissolution of gypsum veins increase the mass equivalent permeability up to 75-300 times depending on aperture width and spacing.

Results and discussion

AJ formation permeability

As mentioned formerly, the weighed average Lugeon value of AJ formation is estimated to be 6 Lu, which is equivalent to about 9 × 10^-7 m/s. Figure 8 shows the variations of this formation permeability against aperture width. Referring to this figure, it is revealed that the dissolution of gypsum veins increase the mass equivalent permeability up to 75-300 times depending on aperture width and spacing. This much increase in permeability advocates the potential hazards of gypsum veins.

Discharge flux inspection

In Figs. 6 and 9, the results of 3D analyses of discharge flux through the dam foundation are shown. In Fig. 6, the dissolution of gypsum is ignored, and the discharge flux is estimated to be 0.075 m³/s. However, in Fig. 9, it is supposed that the gypsum veins finally are washed away, and a passage for seepage flow is evolved. According to this figure, the rate of increase of discharge flux due to dissolution of gypsum veins was obtained to be a third order function of the aperture width. Comparing these figures, the dominant effect of dissolution of gypsum on discharge flux is disclosed. It is seen that as the aperture width reaches to 2 cm due to gypsum dissolution, the discharge flux increases about 240 times. Results of the same analyses with 2D model are shown in Figs. 5 and 11. Referring to Fig. 5 and considering dam length, the initial discharge flux is estimated to be 0.025 m³ / s.

Comparing Fig. 5 reveals that with 2D modeling, the discharge flux increase due to gypsum dissolution is estimated to be around 360 times.

Figure 10 shows the results of variations of discharge flux against the changes of aperture width from 0.04 cm to 2 cm for 3D model. Similar analysis was carried out for the 2D model (typical section), and the results are depicted in Fig. 11. Although the grout curtain is designed to be some 170 m deep, however, complete dissolve of gypsum will severely increase the discharge flux through the foundation ground. As shown in Fig. 10, with 3D model, as the aperture width increases up to 2 cm, the discharge flux rises to 18 m³/s, while with 2D model discharge flux rises to about 0.18 m³/s. For whole dam length, this will be 9.0 m³/s. These figures indicate that the discharge flux with 3D model is about 2.0 folds that of 2D model. This may be attributed to the fact that with the 2D model the lateral components of the seepage velocity vectors are ignored. However, referring to Figs. 5 and 6, this ratio would be about 4 if the gypsum dissolution effect is ignored.

These differences between the results of 2D and 3D approaches may be attributed to the fact that, with 2D seepage analysis, the flow net is planer indeed, and the lateral components of seepage velocity are ignored, so the flow necessarily passes through the successive sections with different permeability. Thus, discharge flux is dominated by the sections with lower permeability. With 3D seepage analysis, however, flow is 3D, and water follows routs with lower energy dissipation, producing higher discharge flux that is realistic.

Pore pressure inspection

In order to inspect the pore water pressure pattern in the embankment and foundation ground, the phreatic water surface through the dam was located using both 2D and 3D models. The results are plotted in Fig. 12. It is seen that the 3D model locates the phreatic surface somewhat higher than the 2D model. This means that the 2D model estimates lower pore water pressure pattern in comparison with the 3D model (Since the total head (H) is accumulation of height head (z) and pore water pressure head (u), so with z being constant in 2D and 3D models, increasing in total head means an increase in pore water pressure). This may be attributed to the fact that with the 2D model, the lateral components of the seepage velocity vectors are ignored. It is noted that in this part of work, the effects of dissolution of gypsum is not taken into account.

Conclusions

In this paper, the seepage behavior of Gotvand Dam considering effects of dissolution of gypsum is evaluated with both 2D and 3D models. Some important conclusions are listed as follows:

1) The 3D analysis estimates the discharge flux about two to four times that of 2D analysis.

2) The 3D model locates the phreatic surface somewhat higher than the 2D model. This means that the 2D model estimates lower pore water pressure pattern in comparison with the 3D model.

3) The dissolution of gypsum veins of AJ formation increase the mass equivalent permeability up to 75-300 times depending on aperture width and spacing; and the rate of increase of discharge flux in Gotvand Dam, due to dissolution of gypsum veins, is a third-order function of the aperture width.

4) According to 3D and 2D approaches, the initial discharge fluxes are 0.075 m³/s and 0.025 m³/s, respectively. However, it is estimated that the gypsum dissolution will increase these figures to about 18 m³/s, 240 times and 9 m³/s, 360 times, respectively.

5) Considering potential hazards of Gypsum veins in foundation ground either deep plastic concrete cut-off wall or a longitudinal gallery for casual grouting in future is recommended.

References

[1] Jeschke A A, Vosbeck K, Dreybrodt W. Surface controlled dissolution rates of gypsum in aqueous solutions exhibit nonlinear dissolution kinetics. Geochimica et Cosmochimica Acta, 2001, 65(1): 27-34

[2] Ford D C, Williams P W. Karst Geomorphology and Hydrology. London: Unwin Hyman Ltd, 1994

[3] Calcano C E, Alzura P R. Problems of dissolution of Gypsum in dam Sits. Bulletin of the Venezuelan Society on Soil Mechanic & Foundation Engineering, <month>July</month>-<month>September</month> 1967, 75-80

[4] Serafim J L, Santos A P, Matos M S. Tridimensional seepage through a dam foundation. In: Proceeding of ICOLD Lausanne. Lausanne: ICOLD, 1985

[5] Uromeihy A, Barzegari G. Evaluation and treatment of seepage problems at Chapar-Abad Dam, Iran. Engineering Geology, 2007, 91(2-4): 219-228

[6] Ghobadi M H, Khanlari G R, Djalaly H. Seepage problems in the right abutment of the Shahid Abbaspour dam, southern Iran. Engineering Geology, 2005, 82(2): 119-126

[7] Ahmed A A. Stochastic analysis of free surface flow through earth dams. Computers and Geotechnics, 2009, 36(7): 1186-1190

[8] Fredlund D G, Rahardjo H. Soil Mechanics for Unsaturated Soils. New York: John Wiley & Sons, 1993

[9] GEO-SLOPE International Ltd. Seep3D User’s Guide. Version 1.1. Alberta: Calgary, 2002

[10] Sadrekarimi J, Saghafi M. Stress dependent permeability effects on ground water flow. Amir Kabir Journal of Science and Technology, 2003, 4(53): 35-42

[11] Freeze R A. Influence of the unsaturated flow domain on seepage through earth dams. Water Resources Research, 1971, 7(4): 929-942

[12] Mahab G C E. Report of International Expert Panel on Upper Gotvand Dam & HEPP. Tehran: Mahab Ghods Consulting Engineers, 2005

[13] Fell R, Mac Gregor P, Stapledon D. Geotechnical Engineering of Embankment Dams. Rotterdam: Balkema A A, 1992

[14] Blyth F G H, de Freitas M H. A Geology for Engineers. 7th Edition. Oxford, UK: Butterworth-Heinemann, 1984, 87

[15] Romanov D, Gabrovsek F, Dreybrodt W. Leakage below dam sites in limestone terrains by enhanced karstification: A modeling approach. Environmental Geology, 2007, 51(5): 775-779

[16] Romanov D, Gabrovsek F, Dreybrodt W. Dam sites in soluble rocks: a model of increasing leakage by dissolutional widening of fractures beneath a dam. Engineering Geology, 2003, 70(1): 17-35

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg