ISSN 00978078, Water Resources, 2012, Vol. 39, No. 7, pp. 811–820. © Pleiades Publishing, Ltd., 2012. Original Russian Text © A.L. Lebedev, A.V. Lekhov, 2011, published in Geoekologiya, 2011, No. 1, pp. 63–74.

RESEARCH METHODOLOGY AND TECHNIQUES

Modeling Changes in the Permeability in a Gypsified FracturedPorous Rock Massif A. L. Lebedev and A. V. Lekhov Geological Faculty, Moscow State University, Moscow, 119992 Russia Email: [email protected] Received November 24, 2009, in final form March 15, 2010

Abstract—A method is proposed for evaluating rock hydraulic conductivity (kx;y, LT–1) in grid model blocks to be used in finitedifference simulation of mass exchange in gypcified fracturedporous rocks containing groundwater. The values of kx and ky were determined taking into account the dominating form of transport (advection or transverse dispersion), contributing most to gypsum dissolution. The method was used to sim ulate groundwater flow in a vertical section of gypsified rocks in the base of the dam at the LowerKafirnigan skii Hydropower System. The changes in rock permeability in the base of the dam were used to evaluate whether there is a need to take seepagecontrol measures (cementgrout curtain, blanket, and a gypsum fill over the blanket). A leaching zone was found to form near the curtain and blanket during the simulation period (100 years) under background seepage conditions. In this period, the total water flow in the lower pool can increase by no more than 10–15%. Under such conditions, the “seepagecontrol” measures can enhance the reliability of the dam at the LowerKafirniganskii Hydropower System. Keywords: gypsum, dissolution, mass transport, dam base, gypsum leaching zone DOI: 10.1134/S0097807812070068

INTRODUCTION In the zone of gypsum occurrence, which is about 7 million km2 in total area [7], changes in the natural conditions, caused by anthropogenic impact on the geological environment, often intensify the mass exchange processes, hence a greater rock permeability and porosity. This problem is especially acute in the hydroengineering construction, where gypsum disso lution and leaching may reduce the stability of rocks. Gypsified rocks are believed to be the least tolerant to the leaching effect [9, 14]. The processes of gypsum dissolution and leaching by groundwater flowing in gypsified rocks (at a scale of a singly fissure, a homogeneous bed or rock massif) have been simulated with the use of advection–disper sion transport equations [1, 2, 9, 13, 14]. Such studies consider the solutions of onedimensional or two dimensional problems obtained by analytical or numerical methods and representing the concentra tion of dissolved gypsum with groundwater flow char acteristics remaining unchanged. The results of grid model simulation of changes in the permeability of a rock massif in a dam base composed of gypsum alone and subject to dissolution are discussed in [16]. Changes in the permeability of a gypsified rock massif as a function of the extent of seepage and leaching pro cesses have received almost no attention of research ers. The objective of this study is to develop method ological approaches to prediction simulation of

changes in the hydraulic characteristics of the model caused by the withdrawal of dissolved gypsum in model blocks with flow assumed twodimensional in the vertical section. The regularities in gypsum occurrence conditions have been mostly determined for areas with karst development. In fractured gypsified rocks, gypsum fills the fracture space between blocks in the form of seams and (or) the pore space within the blocks. Groundwater mostly flows through the fractures that form at the contact between the seams and block walls and, after their dissolution and removal (or when they are absent), in the fractures between blocks. Under such conditions, the schematization of the processes of gypsum dissolution and leaching in model blocks determines the ratios between the mass fluxes of dis solved gypsum in rock blocks and fissures, the rate of mass outflow from blocks into fissures compared with the dissolution rate of seam walls, as well as the length of the path of solution saturation with gypsum. In other words, the time steps and block dimensions are chosen so that the gypsum dissolution in the blocks can be regarded as equilibrium processes. This study gives the results of the final stage of stud ies in a successive investigation of mass exchange between gypsified rocks and groundwater, ranging from a single fracture to a rock massif. The schemati zation of gypsum dissolution and leaching was described in previous works, focused on the kinetics of

811

812

LEBEDEV, LEKHOV

dissolution of a flat surface of natural gypsum in water [5, 6], the rate of gypsum leaching from block pore space [4], and the mass transfer from a single fracture and an element of rock massif at an experimental ground [5]. In the general case, the term leaching refers to the dissolution and transfer into the general body of solu tion of the soluble part of rocks, their porous and insoluble matrix remaining unchanged. In this study, the terms leaching and leaching zone characterize the process of mass exchange between gypsified rocks and groundwater flow at different levels of its studying, including at the scale of rock block, model block, and rock massif. JUSTIFICATION OF THE SIMULATION METHOD The process of changes in hydraulic characteristics caused by gypsum dissolution was simulated using a model of twodimensional flow in a profile, whose thickness is 1 m. The characteristic time of stabiliza tion of flow velocity field is supposed to be much less than the time of changes in the permeability field. Therefore, a steadystate flow problem is solved for the permeability field obtained at each time step. The mass transport equation is written for variable porosity, associated with permeability, in the following form ∂j ∂j ∂ ( nC ) x + y =  + R, (1) ∂x ∂y ∂t where R is mass input due to gypsum dissolution, ML–3T–1; C is dissolved gypsum concentration, C ML–3; ji = –Di ∂ + viC is the mass flux along the ith ∂i coordinate (x or y), ML–2T–1; D is the coefficient of hydrodynamic dispersion, L2T–1; v is the flow rate, derived from the flow problem solution, LT–1; t is the time, T. The flow problem is described by a conven tional twodimensional equation with the transmissiv ity equal to the hydraulic conductivity, multiplied by the unit thickness of the model profile. The mass transport equation is solved by splitting into physical processes, i.e., considering the modeling method, into the process of transformation and trans port [8] ∂j x ∂j y ∂( NC T )  +  =  , (2) ∂x ∂y ∂t ∂( nC R ) – R =  (3) . ∂t In this case, the following conditions are satisfied for the time step Δt = t k + 1 – t k : C T ( t k ) = C ( t k ) transport (initial condition); CR ( tk ) = CT ( tk + 1 ) dissolution (initial condition);

C ( t k + 1 ) = C R ( t k + 1 ) final result. The flow rate under the hydroengineering structure is considerable, so the stabilization time of concentra tion distribution can be supposed to be much less than the time of change of porosity and permeability, the more so as the process of gypsum dissolution proceeds in a narrow zone, because of the high rate of gypsum dissolution. Therefore, the first equation of transport can be assumed stationary, and the realization of mass exchange between water and rock and the recalcula tion of porosity and permeability can be made at the second step of the description of dissolution. The source of mass R is described based on experi mental data on gypsum dissolution rate obtained ear lier [4–6]. Schematization of gypsum dissolution process in model blocks. Until gypsum in the fracture itself is not completely dissolved, its export from the fracture will be minimal because of the small concentration gradi ent for diffusion from the block. Therefore, the pro cess can be conventionally divided into two stages: (i) gypsum dissolution in the fracture, (ii) gypsum export from porous blocks. The changes in model permeabil ity are determined by changes in the permeability of fractures. In this context, the changes in the changes in the total permeability because of dissolution of a part of gypsum in porous rock blocks can be neglected. The solution of the problem of gypsum dissolution in a fracture filled by it, in the case when the flow is possible only along the contact between block walls and the filler, has shown that the transitional zone between the segments of the fracture completely clear and completely filled with gypsum is relatively short [5]. Therefore, we can assume that the passage of gyp sum into groundwater, which is controlled by the undersaturation of the incoming solution and the rate of removal of the dissolved matter, is instantaneous. Thus, before the start of groundwater flow, the model block was a zone of gypsified rocks (zone 2, Fig. 1), which contained gypsum in both fractured and blocks. During groundwater flow, the dissolution of gypsum in fractures is faster than in rock blocks. Therefore, a zone forms in the block, where gypsum is contained only in pore space (zone 1, Fig. 1). Such conditions in the block can be represented as a single zone (1 or 2) or two zones (1 and 2). Such formulation takes into account two most significant transport mechanisms: advection and transverse hydrodynamic dispersion. The change in gypsum mass in an element of the profile (block model) during time dt due to the advection transport along the yaxis is dM C = q y ( C m – C )dt, (4) and that due to dispersion transport along x is ∂C dM D = D T  dt, (5) dy where dMC is the change in the mass of a block, through which unit flow qy passes along the yaxis, WATER RESOURCES

Vol. 39

No. 7

2012

MODELING CHANGES IN THE PERMEABILITY

ML–1; dt is the time interval less than the time of com plete dissolution of the entire gypsum in the element, T; C and Cm are the current and the equilibrium con centrations of dissolved gypsum; dMD is the change of the mass in the element because of transverse hydro dynamic dispersion, ML–1; DT is the coefficient of transverse hydrodynamic dispersion, L2T–1; DT ≈ δTvx; δT is the transverse dispersivity, L; vx is ground water flow rate, LT–1. It is worth mentioning that the profile is assumed to have a unit thickness, which is not taken into account in calculations, as can be seen from the dimension of the unit discharge qy, so the mass has the dimension of ML–1, which, multiplied by the unit length (profile thickness), will yield the required dimension. Clearly, the forecasts of the rate of gypsum dissolu tion in fractures, neglecting its export from blocks, include a safety margin, since the removal of gypsum from rock blocks leads to an increase in its concentra tion in the pore space, thus slowing down its dissolu tion in further downstream in the migration path. Again, this process has a significant effect only in the early time after the dissolution of the entire gypsum in fractures, and its contribution to the total mass flow at the end of gypsum dissolution in the entire block is small. However, to substantiate such assumption requires quantitative estimation of the rate of gypsum export from the block. Schematization of gypsum dissolution process in porous rock block. Let us consider a plate block with a thickness of 2m, L. The process can be assumed long as compared with water migration over the length of the fracture. The increment of the mass in the flow moving along the fracture or the mass discharge of the fracture is equal to the export of mass from the block due to diffusion. Because of the high rate of gypsum dissolution, we assume its concentration on the sur face equilibrium Cm. The equation of mass balance in a fracture free of gypsum at tough conditions of diffu sion in the block is ∂C D ub  =  ( C m – C ), (6) ∂x l where u is the mean water flow velocity in the fracture, LT–1; b is fracture width, L; l is the thickness of the leached part of the block, L; D is diffusion coefficient in the block, L2T–1; C and Cm are the current and equi librium concentrations of dissolved gypsum; x is the migration path, L. Initially, all pores of the block are filled with gyp sum, but dissolution leads to an increase in the thick ness of the zone with open porosity ∂l D ρn b  =  ( C m – C ), (7) ∂t l where ρ is gypsum density, ML–3; nb is block porosity. The product ρnbl is the amount of gypsum exported from a layer with a thickness of l through a unit surface WATER RESOURCES

Vol. 39

No. 7

2012

813

x

MC y

vx

1

Δh

2

MD

Fig. 1. Scheme for evaluating kx and ky in model blocks with Δh taken into account: (1) zone of complete gypsum leaching from fractures (k1, ϕ1); (2) zone of gypsified rocks (k2, ϕ2), MC flow direction at advective transport (ky ⇒ ϕ = ϕ1 + ϕ2); MD is the same at transverse dispersion (kx ⇒ 1/ϕ = 1/ϕ1 + 1/ϕ2).

of the block since the beginning of the process. The increase in the amount of the solid gypsum in the block is equal to the export of dissolvedgypsum mass from the block by diffusion. The boundary conditions C(0, t) = C0 and l(x, 0) = 0 reflect the known concen tration of gypsum in water at the entry of the flow into the fracture and into the block fully filled with gypsum. The solution of this system of equations is 2D ( C – C ) – x  D . t  (8) m 0 ρn b ub From here, it follows that the thickness of the leaching zone linearly decreases along the fracture. In each section x, the rate of its growth gradually decreases. The difference between the thickness of the leached zone in the beginning and the end of the block is constant and equal to xD/(ub), the phase boundary moves down parallel to itself, so that, some time after the beginning of the process, the export of the mass into the fracture in the beginning and end of the block can be assumed practically the same. 1 x C = C 0 +  ρn b D ( C m – C 0 )  . (9) ub 2t The numerical experiment will be made at the fol lowing values of parameters: u = 0.01 cm/s; b = 0.01 cm; D = 1 × 10–6 cm2/s; Cm = 1.7 × 10–3 g/cm3, C0 = 0; x = 20 cm; ρ = 2.3 g/cm3; nb = 0.1. The length of the block lb = 20 cm, m = 10 cm. The values of the thickness of the leached zone in the end of the block and the concentration of dissolved gypsum in the flow moving through a fracture are given l =

814

LEBEDEV, LEKHOV

l, cm; C, g/L 2.0 l

C

1.5 1.0 0.5 0 0.1

1

10 t, year

Fig. 2. Time dependence of l and C at the end of a fracture.

in Fig. 2. The conditions in the block remaining unchanged, all the gypsum will be removed from it within 215 years. After a year, the concentration of cal cium sulfate at the exit from the fracture will be 0.5 g/L; and after 10 years, it will drop to 10% of the concentration of saturation; at the end of the leaching period, it will reach about 2%. The difference between the thicknesses of the leached zone in the end and the beginning of the block is constant and equals to 0.2 cm; after 10 years, it will be less than 10% of the mean thickness. Modeling method. The choice of the modeling pro cedure was based on the following assumptions: the slowest process is the transformation of perme ability distribution in the model space; the concentration field stabilizes faster than the permeability distribution changes significantly; the flow field stabilizes almost instantaneously as compared with the stabilization of the concentration field. Those assumptions allow us to calculate the fields of heads and flow velocities and the field of gypsum concentration in a steady state formulation. After that, the hydraulic conductivity is recalculated and the sim ulation step is repeated. In the recalculation of the permeability, it was assumed that the width of the frac ture does not change after the dissolution of the gyp sum contained in it. The hydraulic conductivity is reevaluated based on the model of fractured medium in the form of orthog onal systems with the same width (b) [12] within a model block. The porosity (n) and the permeability of such medium, for example, in the case of three frac ture systems, can be evaluated as [12] (10) n = 3bFd, 1 3 k P =  b Fd, 6

(11)

where Fd is fracture density, L–1; kP is the permeability (kP = kγ/g , where k is hydraulic conductivity of zone 1, LT–1; γ is kinematic viscosity coefficient, L2T ⎯1; g is gravity acceleration, LT–2), L2. The mass of gypsum in the porous medium of a unit volume was represented as Mn = nρ, ML–3. In the model block (ΔM, ML–1), (12) ΔM = Δx'ΔhM n . The value of Δh was determined from the largest values of MC and MD (i.e., as the most “unfavorable” variant). Substituting into (12) MCorD instead of ΔM, we obtain Δh = MCorD/ ( Δx'M n ). (13) The variables kx and ky were evaluated from the fil tration resistances of the two zones (ϕ1 and ϕ2 in Fig. 1). The obtained values of kx and ky were used at the next step of modeling. The methodological principles proposed in this section for evaluating the permeability characteristics of the blocks of the grid model were used to simulate groundwater flow in a massif of gypsified rocks. The requirements to the natural feature include the solving capacity of the groundwater flow, the occurrence of gypsum in the form of dispersed material of films and seams, filling the fracture space and the adequate knowledge of the region. The region for model appli cation was chosen to be Kafirniganskii anticlinorium, a component of Tadjik depression, where the Lower Kafirniganskii hydrosystem has been constructed since 1984 in the Shaurtauzskii district on the Kafirni gan R., a tributary of the Amu Darya R. The hydrosys tem included a seasonalregulation reservoir 905 mil lion m3 in volume and an HPP with a capacity of 120 thousand kW [3]. The design head differential between the upper and lower pools is 67 m. Now the construction of the hydrosystem is suspended. THE DAM OF THE LOWER KAFIRNIGAN HYDROSYSTEM The base of the dam and its bank abutments are composed of rocks containing many gypsum seams, extending along the channel of the Kafirnigan R. The filling of the reservoir may cause the hazard of intense dissolution of gypsum seams by seeping flow, resulting in a deformation of the hydrosystem. To reduce the rate of gypsum dissolution, it was proposed to use a system of seepagecontrol measures, including grout curtains, asphaltconcrete screen, blanket, and gyp sum filling over the blanket. In the practice of hydroengineering construction, the choice of an optimal variant of seepagecontrol measures is based on their technical and economic comparison (for example, “with a grout” and “with out grout”). Given the criteria of rock seepage stabil ity, the choice of such variant is based on both the groundwater flow characteristics and the cost of con WATER RESOURCES

Vol. 39

No. 7

2012

MODELING CHANGES IN THE PERMEABILITY

struction works. In the absence of criteria, the seep agecontrol measures are used to enhance the reliabil ity of the structure and as a means of additional hydro geological study [10]. No criteria of seepage stability have been developed for the rocks in the base of the Kafirniganskii Hydro system dam. The objective of the use of seepage con trol measures was to protect the gypsum against disso lution and leaching in the dam base (i.e., to reduce the rate of those processes and to prevent chemical suffu sion). Experimental and theoretical studies have been carried out for a unit fracture, filled with gypsum [13]; a homogeneous isotropic porous bed with specified gypsum content [2]; and an experimental ground [11]. In those studies, the need to take seepagecontrol measures was evaluated by the rate and degree of dis solution and/or leaching of gypsum in the dam base [11, 13]. The changes in the hydraulic characteristics of the dam base rocks were also taken into account, though to a lesser degree. The institutes of Soyuzgiprovodkhoz, Tadzhik giprovodkhoz, PNIIIS, VNII VODGEO, and the Geological Faculty of MSU have carried out engi neering geological surveys and works with the aim to justify the seepagecontrol measures in the construc tion area. Based on the results of those studies, we give below a concise characteristic of the geological and hydrogeological conditions in the work area and the flow scheme of the rock in dam base. A BRIEF CHARACTERISTIC OF THE GEOLOGICAL AND HYDROGEOLOGICAL CONDITIONS IN THE CONSTRUCTION REGION OF THE DAM OF THE LOWERKAFIRNIGAN HYDROSYSTEM In the construction area, the river has a wide, well developed valley with stepwise profile, filled with a stratum of alluvial deposits of kafirniganskaya suit of Miocene, which compose river floodplain and for abovefloodplain terraces. Depending on the litholog ical features, the deposits in the area can be divided into three units: the lower (N1kf1), the middle, (N1kf2), and the upper (N1kf3). Lithologically, these are inter calating beds of reddishbrown, brown, and gray, polymictic, coarse and mediumgrain sandstones, as well as argillites and aleurolites with the predominance of pore and film types of cement with gypsum, carbon ate–gypsum, and iron–clay composition. Tectonically, the territory is an asymmetric syn cline, extended from north to south, with counter inclination of rocks under the angles of 45–55°. Three systems of fractures, coinciding with the direction of faults are predominant (the first one, along the rock bedding). The rock massif has a heterogenic–block structure with block sizes varying from 18 to 32 cm in fault zones and from 30 to 62 cm beyond those WATER RESOURCES

Vol. 39

No. 7

2012

815

zones (i.e., δ1 ≈ 0.2–0.5, δ1 is the longitudinal disper sivity, L). Gypsum in the rocks occurs in two generations: (I) as a filler of fractures, monomineral; (II) dispersed and superficial, developed in the form of film and con tact cement, inclusions, pockets, etc. The upper part of the section of primary deposits (down to 40 m) is intensely weathered, leached, and washed. The frac tures below this zone are filled with white crystalline gypsum. Two types of groundwater are identified: uncon fined pore water of Quaternary deposits and fracture– pore waters of bedrocks, which form a single aquifer. Beyond the zone of modern impact of exogenous pro cesses (the depth greater than 8–25 m), the kafirni ganskie deposits feature a relatively low background permeability—the specific discharge is less than 0.003 L/(min m2), which is nearly zero. In the top part of the section, the specific discharge increases 2– 3 times up to 0.014–0.017 L/(min m2), i.e., the rocks can be regarded as lowpermeability. The zones of higher permeability occur locally and are generally associated with tectonic faults, higher fracture density, and gypsification. The difference between the back ground and local permeability is estimated at 2– 3 orders of magnitude. Kafirnigan R. water in its lower reaches shows higher mineralization (0.5 g/l). Subsoil water mineral ization is 1.8–20 g/L (pH 7.3–7.7; 20°C), and its chemistry corresponds mostly to sulfate–sodium sub type of sulfate water type (Kurnakov–Valyashko clas sification). The results of calculation of ionic equilib ria have shown that (i) subsoil water is undersaturated with respect to gypsum to a different degree (0.5– 1.8 g/L); (ii) during the mixing of subsoil water with Kafirnigan R. water, the gypsumsolving capacity of the solution increases, reaching its maximal value (Cm = 3.87 CaSO4 ⋅ 2H2O g/L) in “pure river water”, i.e., the replacement of the fracture water by water from the reservoir will not lead to the formation of solutions more corrosive with respect to gypsum. According to the results of previous studies [4–6, 15], the process of mass withdrawal from fractured gypsified rocks can be regarded as a series of successive processes of mass withdrawal from fractures filled with gypsum, and leaching the latter from the pore volume of rock blocks. The model of fractured medium repre sents it as three orthogonal systems with the same frac ture widths b within a model block. With the minimal linear dimensions of the model blocks (Δx', Δy') of 5 m, the value of ΔM can be calculated taking into account the equilibrium characteristics alone. With the gypsum content of the rock in the block of ~5%, the values of diffusion coefficient D of Ca2+ ions vary within (3.4–7.7) × 10–6 cm2/s (25°C, bidistilled water) and (0.5–1) × 10–6 cm2/s (20°C, Kafirnigan R. water) at gypsum concentration not higher than 1–2% [4–6, 15].

816

LEBEDEV, LEKHOV C D E

305 m k1

1 2 2

3

4

5

3 777

4

5 1.33

50 m

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

385 m B

A

25 m

330 m

k2

6

6 7 8 9 1011121314151617 18192021222324252627 28 29 30 31 32 33

34

35

36

Fig. 3. Schematic model grid and streamlines in modeling groundwater flow in dam basement of the NizhneKafirniganskii Hydraulic System: (1) leaching zone; (2) grout curtain with hydraulic conductivity of k3; (3) gypsum filling; (4) impermeable boundary; (5) streamline at Δt = 0 without curtain (constructed with the help of discharge profiles); (6) block numbers. BCDE contour is dam body. AB is the blanket.

MODELING CHANGES IN THE PERMEABILITY OF DAM BASE CAUSED BY GYPSUM LEACHING The need for the use of seepagecontrol measures was assessed by the results of simulation of water flow in dam base with different arrangement and length of grout curtain and gypsum filling over the blanket. A scheme of twobed structure with an impermeable top layer was used (see Fig. 3). Changes in the hydraulic characteristics of the medium were regarded as the result of gypsum leaching alone, so the issues of gen eral ecological expertise were not considered, as well as physicochemical and mechanical suffusion. The results of calculations showed that during 50– 100 years of dam operation, gypsum can be completely leached from rock blocks (D = 1 × 10–6 cm2/s, nb = 0.02–0.05). It can be supposed that an increase in porosity due to the dissolution of 2–5% of gypsum will not have significant effect on the increase in the per meability of dam base, since the permeability of rocks is mostly determined by the permeability of fractions, so at this stage of modeling, only the leaching of gyp sum from fractures was evaluated. The characteristics of the object were chosen to be the dependence of the total groundwater flow (Q, L3T–1) and the total mass output of dissolved gypsum (M, M) in the lower pool on the simulation time. The domain of the model with fractures containing gypsum was spec ified by a boundary condition of concentration Cm. After the complete dissolution of gypsum in the block, this boundary condition was withdrawn. A boundary condition on the initial concentration was specified on the reservoir bed in accordance with river water chem istry (C) or the presence of gypsum filling (Cm). The boundary condition on concentration was not speci fied in the lower pool. The variables Q and M were evaluated by the following procedure. 1. Source data inputting. At this stage, the grid model blocks were assigned permeability characteris

tics (k1 = 0.25–1 m/day; k2 = 0.01–0.05 m/day; k3 = 0.01 m/day), the concentrations of dissolved gypsum (Cm = 3.87 g/L; C = 0.16 g/L) and the dispersivity (δ1 = 0.2–0.5, δT = 0.1δ1). Gypsum filling and grout curtains were specified by the values of Cm in the blocks N (8–9)/1 and k3 in blocks N 18/(1–5); N 10/(1–3) (the top numbers are the numbers of blocks in the hor izontal direction, the bottom number is that in the ver tical direction). The size of the blocks Δy' = 5 m (N 1– 15); Δx' = 50 m (N 1–5, 34–36), 35 m (N 33), 30 m (N 29–32), 25 m (N 28), 20 m (N 6–27). 2. Calculating the fields of head (H) and concentra tion (C) was made by using the direct Gauss–Jordan method with the help of TRANSFER program. 3. Evaluating Q and M. The value of Q was calcu lated from Darcy law at the interface between blocks N 1/(29–36) and N 2/(29–36) (see Fig. 3). The values of Qb (L3T–1) for each block thus obtained were summed and the total value of Q was used to plot the dependence of the form Q = Q(t). The values of (M) were calculated based on the difference between the values in the lower (Mb, M) and upper (Mu, M) pools. Ml,u = QbΔtC, where C is the concentration in the ith block. 4. The length of the zone of complete dissolution of gypsum (Δh) was calculated by equations 4–5 and 10– 13. The values of Fd were calculated from the fractur ing modulus 1–7 [3, 11], the results ranging from 0.01 to 0.09 cm–1. At k1 = 0.25 m/day and considering the most frequent values of Fd (0.028–0.85, we have b = (3.7–2.5) × 10–3 cm and n = 0.0003–0.0007. In the calculations, we used the largest value n = 0.0007 (k1 = 0.25 m/day). The values of n at other values of k1 were calculated similarly. 5. Evaluating kx, ky. The characteristics specified for the given model block in the grid model included the initial values of Δt, n, Δx', Δy', the differentials of the heads and hydraulic conductivities between the blocks. The results of calculations (kx or ky) were WATER RESOURCES

Vol. 39

No. 7

2012

MODELING CHANGES IN THE PERMEABILITY Q, m3/day 0.96 1 3

2

0.92

0.88

0.84 0

40

20

80

60

100 t, year

Fig. 4. Dependence of Q on t and Δt without curtain (k1 = 0.25 m/day, k2 = 0.01 m/day, n = 0.007 and δ = 0.2): (1) Δt = 40 years; (2) Δt = 10 years; (3) Δt = 5 years (within the range from 0 to 50 years) and 7 years (at 50–100).

assigned to the ith block of the grid model in the leach ing domain (see Fig. 3). The length of the grout curtain (LZ) was evaluated taking into account changes in the total flow in the lower pool (ΔQ , L3 T–1) with the grout curtain and without it (Δt = 0; k2 = 0.01–0.1 m/day). The maximal

7

8

9

10

11

12

25

value of ΔQ (7%) was obtained at LZ = 25 m (k2 = 0.01 m/day). The value of ΔQ showed practically no changes at farther increase in LZ. So, the use of a grout curtain will be most effective when its length is equal to the thickness of the first, most permeable, layer, as is commonly accepted in the choice and designing of seepagecontrol measures. The time step (Δt) was chosen at the stage of pre liminary calculations “without the grout curtain” both at fixed values of (5, 10, 20 years), and at variable val ues (from 2 to 30 years). At the runs of different vari ants, the values of Q stabilized at Δt ≤ 5, so the maxi mum admissible value of Δt was chosen to be 5 years for the first half of the calculation period, and 7 years for the second half (Fig. 4). MODELING RESULTS The process of gypsum leaching is most intense in the zones of the roof of the gypsified rock layer, within which the transport is mostly due to advection (i.e., in the zone where streamlines cross the interface between groundwater unsaturated and saturated with gypsum at nearly right angle). In such areas (blocks N (6– 7)/(7–9, 16–18), Fig. 5), specific holes will form, i.e., the washing out of gypsum is mostly directed down

Block numbers 13 14 15

16

17

18

19

20

21

1

26

2 3

27 Gypsum occurrence depth in fractures, m

817

28 4

29

5

30 31

6

32

7

25

1

26

2

27

3

28 29

4

30

5 6

31 32

7

Fig. 5. Variations in gypsum occurrence surface in fractures during water motion in the rocks of dam basement: over (1) 10 years, (2) 20, (3) 30, (4) 40, (5) 60, (6) 80, and (7) 100 years. The top and bottom plots are for the cases without and with the curtain. The dashed line is for advection transport, the full line is for dispersion transport. WATER RESOURCES

Vol. 39

No. 7

2012

818

LEBEDEV, LEKHOV Q, m3/day 4

Q, m3/day 0.97 1

0.93

3

1

0.89 2 2

3

0.85 0.81

1 3

2

0.77 0

0

20

40

60

80

40

20

60

80

100 t, year

100 t, year

Fig. 6. Dependence of Q on t and k1 without curtain (k2 = 0.01 m/day; δ = 0.2); (1) k1 = 1 m/day, n = 0.001; (2) k1 = 0.5 m/day, n = 0.0008; (3) k1 = 0.25 m/day, n = 0.007.

Fig. 7. Time dependence of Q at k1 = 0.25 m/day, k2 = 0.01 m/day, n = 0.0007, and δ = 0.2; (1) without curtain; (2) with a curtain; (3) with filling.

ward. The maximal increment in flow (Q) takes place in the first 5–10 years (Figs. 4, 6–8), i.e., when leach ing takes place on the surface of the roof of gypsified rock layer. Within this period, the use of one curtain reduces Q by ~5–7% (Figs. 7, 8 a, b), while two cur tains reduce it by ~11%. The results of modeling have shown that at leach ing processes of such character, the value of Q during the simulation period without seepagecontrol mea sures increased by only 10–15%. The use of the cur tain (LZ = 25 m) became inefficient as early as the first half of the simulation period (see Fig. 7), though it increased the total amount of dissolved gypsum (M). Overall, in the simulation period, 47000 and 60000 t of gypsum can be leached from dam base rocks (at the section length of 1500 m) without and with the cur tain, respectively (see Fig. 9).

domain in the upper pool) coincides with the direction of streamlines (see Figs. 3, 5). The similar character of the processes of gypsum dissolution and leaching in the rocks composing the dam of the NizhneKafirni ganskii Hydrosystem was demonstrated by B.S. Sherzhukov, A.S. Malyshev, and N.K. Golog vanov [13]. In their study, they considered a one dimensional problem of advection diffusion in a single Q, m3/day 1.18

(a) 1

1.16 1.14

2

1.12

DISCUSSION OF RESULTS The values of Q and M obtained at the use of differ ent seepagecontrol measures correspond to certain flow conditions (k1 = 0.25 m/day; k2 = 0.01– 0.05 m/day; n = 0.0007–0.003; δ1 = 0.2–0.5), which can be assumed background for dam base rocks and at which, the use of “protection means” (including gyp sum fill) is not reasonable, since the effect of the leaching on the increase in water flow in the lower pool is quite insignificant (~10–15%). In higher permeability zones (k1 > 0.25 m/day, k2 > 0.01 m/day) the values of Q and M can be large. In such cases, the application of seepagecontrol mea sures should be based on the size of zones, their hydraulic characteristics, hydrosystem operation time, and the locations and sizes of the curtain, blan ket, and fill. The results of modeling have shown that the domi nating direction for the extension of the domain of gypsum leaching in dam base (hereafter, the leaching

1.10 1.08 (b)

0.92

1

0.88

2

0.84

3

4

0.80 5

0.76 0.72 0

20

40

60

80

100 t, year

Fig. 8. Time dependence of Q (a) at k1 = 0.25 m/day, k2 = 0.05 m/day, n = 0.003, and δ = 0.2: (1) without curtain; (2) with a curtain; (b) at k1 = 0.25 m/day, k2 = 0.01 m/day: (1) without curtain, n = 0.003, δ = 0.2; (2) with filling, n = 0.003, δ= 0.2; (3) with a curtain, n = 0.004, δ = 0.5; (4) with a curtain, n = 0.003, δ = 0.2; (5) with two curtains, n = 0.003, δ = 0.2. WATER RESOURCES

Vol. 39

No. 7

2012

MODELING CHANGES IN THE PERMEABILITY

fracture filled with gypsum. The obtained solution in the form of the concentration distribution of dissolved gypsum for a model section in the dam was used to predict the dissolution of seams with a total thickness of (2.86 m). It was found that, at the water undersatu ration with gypsum by 1.5 g/L, a hole up to 1 m in depth can form in the upper pool; this fact is quite comparable with the results of this study (i.e., at the undersaturation of 3 g/L, a hole with a depth of ~4– 5 m will form, see Fig. 5). In [16], a grid model was used to simulate the effect of dissolution of fracture walls on changes in the per meability of a rock massif completely composed of gypsum. A profile problem with water flow in dam base with a grout curtain was solved. The model domain (750 × 375 × 1 m) was divided into impermeable blocks (7.5 × 7.5 × 1 m) by two systems of mutually perpen dicular fractures with the same initial value b = 0.01 cm. The rate of gypsum dissolution was chosen for each fracture depending on Q and b. The values obtained with the use of equations of diffusion kinet ics, heterogeneous chemical reactions on gypsum reaction surface were compared, and the lesser values were taken. The model results were used to construct the dependence of Q on t for the lower pool [16]. It was found that, as the gypsum dissolution domain extends in the direction from the upper to lower pool, the value of Q increases by ~7–15% (i.e., similarly to our calcu lations) with b values up to ~2–4 cm. When the boundary of this domain reaches the lower pool, the value of Q abruptly increases by several orders of mag nitude at b of up to 20 cm. The results of these studies and the studies [13, 16] show that the dissolution and leaching of gypsum near the upper pool have no significant effect on rock per meability in dam base rocks. However, the extension of the leaching domain to the lower pool can cause an increase in Q by several orders of magnitude. Such character of changes in flow rate was observed in theoretical and experimental studies of mass trans fer from a single fracture with walls composed of gyp sum. Thus, when the fracture is longer than the path of water saturation with gypsum and with the values of b and Q typical of rock massifs, the dissolution takes place near the entry and the value of Q shows almost no changes. However, the extension of the dissolution domain from the entry toward the exit section of the fracture leads to a considerable increase in the flow rate [9, 14, 16]. Thus, in the justification of seepagecontrol mea sures for gypsified fracture–porous rocks in the base ments of hydrostructures, it is reasonable to consider the possibility of driving the flow via the longest path in the domain of the upper pool. Possible sources of errors and the perspectives of their taking into account in the future studies. The schematization of gypsum leaching in model blocks WATER RESOURCES

Vol. 39

No. 7

2012

819

M, mg 6.0E + 07

2

4.0E + 07 1 2.0E + 07

3 0.0E + 00 0

20

40

60

80

100 t, year

Fig. 9. Time dependence of the total mass of gypsum dis solved in dam basement (M) at k1 = 0.25 m/day, k2 = 0.01 m/day, n = 0.007 and δ = 0.2: (1) with a curtain, (2) without curtain, (3) with filling.

assumed its instantaneous dissolution with the forma tion of only two zones: complete leaching and gypsi fied rocks (see Fig. 1). The evaluation of kx and ky was carried out for the conditions corresponding to the leaching of the maximal amount of gypsum, i.e., the appropriate transport mechanism (advection or trans verse dispersion) was chosen. Such simplifications with the aim to use the most unfavorable parameters are admissible when running different variants of protection measures, since they allow forecasting the changes in rock permeability in dam base with a sure safety margin. A more accurate estimate of changes in the hydraulic characteristics will require more detail schematization of the pro cesses of gypsum dissolution and leaching in transport mechanisms. The results of simulation over the model period (100 years) showed the process of gypsum leaching from dam base rocks to proceed only in a zone at the upper pool (see Fig. 3). Gypsum dissolution in a zone near the lower pool can be accompanied by an abrupt increase in Q by several orders of magnitude, as shown in [16]. Therefore, when calculating groundwater flow in such studies, one should not limit the calculation to a specified time period, but carry out this, taking into account the processes of gypsum dissolution and leaching all over the flow path. CONCLUSIONS In the finitedifference modeling of mass exchange processes in massifs of gypsified fractured and frac tured–porous rocks containing groundwater, the cal culations of model block permeability can be based on the mechanism of the predominant transport (advec tion or transverse hydrodynamic dispersion), which accounts for the largest amount of dissolved gypsum. This creates a safety margin in problem solution.

820

LEBEDEV, LEKHOV

When groundwater not saturated with respect to gypsum flows through a horizontal bed of gypsified fractured–porous rocks, two zones of successive phases of gypsum leaching can be identified: on the surface of the roof and within the layer. In the former zone, the transport is due to advection and transverse dispersion, while in the latter zone, it is mostly due to advection. During gypsum leaching in the zone near the upper pool of the dam, an increase in the permeability of the massif of gypsified fracture–porous rocks can be very small. Under such conditions, seepagecontrol mea sures are a means for increasing the reliability of the hydraulic structure (rather than a means to prevent seepage losses), as was the case with the seepage in the base of the anticipated dam of the NizhneKafirnigan ski hydrosystem. REFERENCES 1. Verigin, N.N. and Sherzhukov, B.S., Diffusion and Mass Exchange during Liquid Flow in Porous Media, in Razvitie issledovanii po teorii fil’tratsii v SSSR (1917⎯1967) (Progress in the Theory of Fluid Flow in Porous Media in the USSR (1917–1967)), Moscow: Nauka, 1969, pp. 237–313. 2. Dem’yanova, E.A., Gypsum Dissolution and Export by Water Flowing through Dam Basements, Inzhenernaya Geologiya, 1986, no. 6, pp. 23–33. 3. Kornev, Yu.P., Kondrat’ev, A.F., and Gavrikov, B.I., NizhneKafirniganskii Hydrosystem, Gidrotekh. Melior., 1985, no. 7, pp. 11–17. 4. Lebedev, A.L., Lekhov, A.V., Sokolov, V.N., and Svi toch, N.A., Gypsum Leaching Rate from Sandstone Pore Space,Geoekologiya, 2003, no. 5, pp. 438–447. 5. Lebedev, A.L., and Lekhov, A.V., Mass Export from Gypsified Fractured Rock with Groundwater, Water Resour., 1999, vol. 26, no. 3, pp. 277–285. 6. Lebedev, A.L. and Lekhov, A.V., Dissolution Kinetics of Natural Gypsum in Water at 5–25°C, Geokhimiya, 1989, no. 6, pp. 865–874. 7. Maksimovich, G.A., Osnovy karstovedeniya (Funda mentals of Karst Science), Perm: Perm. kn. izd., 1963, vol. 1.

8. Mironenko, V.A. and Rumynin, V.G., Problemy gidro geoekologii. Tom 1. Teoreticheskoe izuchenie i modeliro vanie geomigratsionnykh protsessov (Problems of Hydroecology. Vol. 1. Theoretical Studying and Model ing Geomigration Processes), Moscow: Izd. MGTU, 1998. 9. Oradovskaya, A.E., On the Methods of Studying the WaterTolerance of SemiHard Rocks GypsumBear ing Rocks, in Problemy inzhenernoi geologii v stroi tel’stve. Mater. soveshch. po inzhenernogeologicheskim issledovaniyam skal’nykh osnovanii gidrotekhnicheskikh sooruzhenii (Problems of Geoengineering in Construc tion. Mater. Conf. on Geoeng. Studies of HardRock Basements of Hydraulic Structures), Moscow: Gos stroiizdat, 1961, pp. 182–195. 10. Pavlovskaya, L.N., and Fomina, N.E., Studying the Effect of the Elements of Subsurface Contour of a Con crete Dam on Counterpressure, Izv. VNIIG, 1981, vol. 146, pp. 46–50. 11. Pokrovskii, G.I., Kol’tsov, G.A., Golovko, A.N., and Tregubko, V.P., Safety Assurance of the Dam at the NizhneKafirniganskii Hydrosystem, Melioratsiya i Vodnoe Khozyaistvo, 1988, no. 9, pp. 14–16. 12. Romm, E.S., Fil’tratsionnye svoistva treshchinovatykh gornykh porod (Hydraulic Properties of Fractured Rocks), Moscow: Nedra, 1966. 13. Sherzhukov, B.S., Malyshev, A.S., and Golovanova, N.K., Forecasting Bed Gypsum Dissolution in the Basements of Hydraulic Structures, in Metody rascheta protsessov massoperenosa v gidrogeologicheskikh issledovaniyakh (Methods for Calculating Mass Transport Processes in Hydrogeological Studies), Moscow: VNII VODGEO, 1984, pp. 26–28. 14. James, A.N., and Lupton, A.R., Gypsum and Anhy drite in Foundations of Hydraulic Structures, Geotech nique, 1978, vol. 28, no. 3, pp. 249–272. 15. Lebedev, A.L., Gypsum Dissolution and Leaching Pro cesses: The Influence on Permeability and Porosity of Fractured–Porous Rock Massif, Proc. Int. Symposium. Engineering Geology and the Environment, Athens, Greece, 23–27 June, 1997, Balkema, pp. 823–827. 16. Romanov, D., Gabrovsek, F., and Dreybrodt, W., Dam Sites in Soluble Rocks: A Model of Increasing Leakage by Dissolutional Widening of Fractures Beneath a Dam, Engineering Geology, 2003, vol. 70, p. 1735.

WATER RESOURCES

Vol. 39

No. 7

2012