- Email: [email protected]
Particles’ transport in a single fracture under variable flow regimes
ADES 431
Advances in Engineering Software 30 (1999) 327–337
Particles’ transport in a single fracture under variable flow regimes C. Masciopinto* National Research Council, Water Research Institute, Via F. De Blasio 5, 70123 Bari, Italy Received 1 June 1996; received in revised form 1 October 1998; accepted 1 October 1998
Abstract A mathematical model to account for non-laminar terms, during simulated transports in fractured rocks has been developed. The relevant conceptual model consists of a two dimensional grid block equivalent to a variable fracture plane. The block aperture variability is expressed by means of stochastic procedures, based upon some characteristics of fractured media. One of them, e.g. the equivalent hydraulic aperture, may be estimated by means of the results of field tests. The solution of Navier–Stokes and Darcy–Weisbach’s equation for steady state flow between two parallel plates with constant aperture allows the systems of equations under respetively, laminar and non-laminar conditions, to be wirtten. The flow model is a standard finite-difference method with a modified version of Darcy’s equation to account for non-laminar terms that are important at high water velocity. To simulate solute transports, the particle tracking technique was applied. 䉷 1999 Elsevier Science Ltd. All rights reserved. Keywords: Fractured media; Radial flow; Non-laminar regime; Numerical model; Fracture apertures
1. Introduction Particle tracking technique can be very useful for the identification of the control area in well field. In this article such a technique is used to study the tracer transport in a square area of a fractured aquifer. The study of pollutant transport in an aquifer cannot neglect the physical characteristics of soils, being either fractured or granular. The usual diffusive model, which describes flow and transport of pollutants in terms of an ‘‘equivalent porous media’’, may be inconsistent in fractured systems. In fact, in systems with few fractures or in a heterogeneous porous media, the basic assumption of the continuum approach, i.e. the Representative Elementary Volume (REV) to which all physical properties of the medium are referred, is often missing. Recently the efficacy of the stochastic approach to study heterogeneous porous formations has been showed [4,35]. Medium heterogeneity can be described in a schematic representation, also three-dimensional, by means of significant parameters [12]. In rock systems with few fractures some authors [8,21] proposed fractal theory, to reproduce a single fracture and, to study the fluid flow; the particle tracking can be applied, as clearly explained by Schwartz et al. [36], to study pollutants transport in 2- or 3-D scale [1,29,42,44]. This method * Corresponding author. Fax: 0039 80 5313365. E-mail address: [email protected] (C. Masciopinto)
was applied by Schwartz et al. [36] to study solute transport in fracture sets having constant apertures and different orientations. For Neretnieks [31] and Tsang and Tsang [40], a single fracture can be considered equivalent to a set of independent channels, with variable cross-section. The apertures can be derived by log-normal distributions [38] of values and can be generated by means of an exponential function of the spatial covariance of the apertures. Fracture apertures can be obtained when the values of mean, standard deviation and correlation length are known. Moreno et al. [27] applied this technique to study the dispersion through a single fracture with variable aperture under uniform and laminar one-dimensional flow conditions. A computational code has been developed here, to model non-uniform steady state flow under variable flow regimes in a variable-aperture fracture plane (Fig. 1). These groundwater fluid flow conditions are frequently verified inside the capture zones of pumping wells. The flow model is based on a standard finite difference method with a modified version of Darcy’s equation, to account for non-linear terms that are important at high velocity values. An application of this model allows the evaluation of the non-laminar effects on the pollutant transport in fractures to be done. Clogging effects and related phenomena have not been considered in this article. The fact that the fracture aperture is as low as 1–20 mm, in the narrowest parts of many fracture channels, means that most of the soil particles that can be moved, belong to the clay fraction. Clay minerals tend to
0965-9978/99/$ - see front matter 䉷 1999 Elsevier Science Ltd. All rights reserved. PII: S0965-997 8(98)00092-1
328
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Fig. 1. Fracture generation parameters for simplified models.
be stuck by interacting with fracture minerals. Thus, the drug forces need to overcome bonds between transported soil particles and fracture minerals. To avoid clogging effects caused by transported soil particles a flow rate exceeding a minimum value (about 8.6 m/d [10]) is required. These effects will be considered in the next article.
2. Fluid flow motion Fluid flow can be studied applying the solution of the Navier–Stokes’ equation for average velocity in a single fracture of constant width, b, bounded by two parallel impervious plates [39]
g 2 b 7f; V ⫺ 12m
1
where m and g are respectively the fluid dynamic viscosity and specific weight (g /m is about 10 7 m ⫺1s ⫺1 for water at
20⬚C), while f is the piezometric head in the extremities of the channel with constant length Dx and constant crosssection Dy × b. Using a finite difference method with a block centred approach, the grid nodes (Fig. 2) may be defined using a grid discretization step size capable to replace flow streamlines with subsequent segments in x and y direction. Moreover, if i and j are two consecutive grid nodes, oriented channels in x or y direction can be defined. As the flow is parallel to the grid, the effects of grid orientation on the flow velocity [19] are not significant [28]. These channels can be considered with constant length Dx , or Dy, and constant cross-section Dy × b, or Dx × b. Discharge values can be obtained from the equation Qij
g b3 Df Dy Dx m 12
2
Dx and Dy were selected short enough so that one-dimensional flow over each channel occurred and the relation
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
329
Fig. 2. Grid nodes and channels.
between interconnected piezometric heads and discharge is given by
fi ⫺ fj
Dfi ⫹ Dfj Qij 2 1=12 g=mb3i Dy 2=Dx Qij 1=12 g=mb3j Dy 2=Dx !# 1 1 ⫹ 3 Qij Rij ; b3i bj
Qij 6
m Dx g Dy
X
Qij
J
⫹
"
grid, without sink and source, we can write
3
where Rij is the flow resistance of the channel between nodes i and j under laminar regime. Eq. (3) proposed a linear gradient interpolation between adjacent blocks with constant apertures bi and bj. In fact, using bi and bj into Eq. (2), two piezometric head differences can be found, Df i and Df j, respectively. Thus, because of nodal step changes in thickness, the linear interpolation of the gradient gives a channel conductance, 1/Rij, equal to the harmonic mean of the nodal values. To preserve the net discharge across each node of the
X fi ⫺ fj 0; Rij J
4
where the converging channels, J, are three, instead of four, along the borders of the studied area. By using Eq. (4) for all grid nodes a system of linear equations can be obtained; the system solution can be found after imposing suitable boundary conditions. To improve this method, the flow resistances were corrected to consider inertial forces for high Reynolds numbers (⬎ 10). Under these field conditions, the Fanning factor [3] f can be used to evaluate the frictional losses inside each channel f
2g Df ; V 2 Dx
5
where g is the acceleration of gravity. The Fanning friction factor in laminar condition, for fluid flow motion between two parallel impervious surfaces, can be derived combining
330
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Fig. 3. Relationships between the Fanning friction factor and Reynolds number.
Eqs. (1) and (5). If the hydraulic radius Ri is Ri
Dy × b b 艑 ; 2Dy ⫹ 2b 2
6
valid for b ⬍⬍ Dy: Fanning values can be obtained as f 24=Re, where Re
4rVRi b 2rV m m
7
is the Reynolds number for channel with constant rectangular cross-section [24]. Also for laminar flow through porous (granular) media, the relationship between Fanning friction factor and Reynolds number gives a straight line portion, in a log–log plane, expressed by f 1000=Re. Thus under laminar conditions the f f Refunctions, are two parallel straight lines in a log–log plane for both porous media and impervious surfaces. Under non-laminar conditions, in agreement with Gale [17], measurements of the friction factors are rare, despite the importance of such values in hydrological studies. Results obtained by Blumberg and Curl [7] and Benedini et al. [5], are not suitable in case of limestone that does not have very large conduits (caves) or regular channels (e.g., bricks). Assuming that the shape of relations f f Re holds also for non-laminar regimes, the Fanning friction factor for flow through a single fracture can be obtained by a translation of the granular porous media data, using the constant 24/1000 (Fig. 3). In this way the Darcy–Weisbach friction factor, l 4·f ; can be evaluated by means of the best-fit of porous media data [34], shifted by value 96/1000 105:12 992:31 300:70 96 ⫺ : 8 l 15:53 ⫹ p ⫹ Re 1000 eRe Re This equation gives a best fit also for laminar (Re ⬍ 2) and
turbulent [22] (Re ⬎ 10000) regimes. To verify this behaviour, a forced tracer test with high discharges was carried out and some points under non-laminar regime were found [16] (see Fig. 3). Using the Darcy–Weisbach equation, the corrected discharge can be obtained by s 4bg Df Qij b × Dy 9 × l Dx after that, the solution of the system (4) was found and the Re numbers calculated by means of the velocity values given by Eq. (1). The relation between interconnected piezometric heads and discharge, under generic flow regime, is now given by ! Dfi ⫹ Dfj Q2ij li;j 1 1 × fi ⫺ fj ⫹ 3 2 8gDy Dy=Dx b3i bj " !# l i;j Dx 1 1 × ⫹ 3 Q ij 2 D y b3i 8g D y bj Q ij
2
R 0ij ; 10
Rij0
is the corrected flow resistance of the channel where between nodes i and j. Without sinks and sources, in each grid node, we can write s X fi ⫺ fj X 0 11 Qij R 0ij that represents a system of non-linear equations, the solution of which can be found after imposing suitable boundary conditions. The developed computational code employs the successive over relaxation method (SOR) to solve
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
systems (4) and (11), after linearization. When laminar conditions occur, systems (4) and (11) give the same solution.
3. Tracer transport The phenomenon of solute transport is investigated by tracking the particles advected through an ideal single fracture with variable apertures and with a behaviour which is equivalent to that of a fractured media (see Fig. 1). The particles can be introduced into a generic node of the domain. After that, the constant discharges within each channel were evaluated. However, during pumping simulations, the particles have been left in border nodes and collected in the central node of the domain. Particles coming to an intersection (i.e., the generic block-centred grid node) are distributed in the outlet branches with a probability proportional to the flow rate [2,33]. This method is consistent with an instantaneous and complete solute mixing in each node and the majority of particles take the flow path with higher discharge. If dM nC0 x0 dx) is the polluter mass with concentration C0 at t t0 at the inlet of the channel, its Lagrangian coordinates can be expressed as a function of the probability density of the flow path particles [13,32]. Here, the numerical solution was used. Following each particle along its pathway as the Lagrangian point of view requires, it migrates inside channels, oriented in the x or y direction (see Fig. 2), consecutively. Within each channel uniform and steady-state flow occurs and for preserving the discharge, as the channel cross section is constant, also the velocity must be constant (incompressible fluid and non-deformable rock). To observe these hypotheses the size of the grid block cannot be very large. As the mass flow at outlet location is proportional to the number of particles that reach that location, the concentration is derived from the number of particles per unit time that arrive at a collection point divided by the flowrate at the same location [29]. The mean residence time of the particle inside each grid element of constant length Dx and constant cross-section Dy × bij (i.e., channels in x direction) can be determined by using Dx ti ; Vij
12
where Vij is given by arithmetic mean of the apertures (see Fig. 2). As Dx Dy the same equation can be applied for channels oriented in the y direction. This equation was preferred to that proposed by Moreno et al. [27] to avoid lateral effects on the volumetric flux from node i to node j [20]. The total particle residence time is given by the incremental addition of these times over the entire particle path from the inlet to the outlet position. Following Moreno et al. [27], under equilibrium conditions, a sorbing tracer in each channel can be considered
331
using the retardation factor, Ra, for each node Ra 1 ⫹ 2
Ka ; bi
13
where Ka is the surface distribution coefficient which can be obtained from static sorption experiments (i.e., batch tests) [9,43]. In fractured systems with large macropores, the tracer diffusion into fracture walls cannot be neglected and kinetic equations should be adopted to consider the effects of diffusion. 4. The fracture generation The computational code has been written to solve fluid flow and pollutant transport in the plane of the statistically generated fracture. To solve the fluid flow problem, steady state boundary conditions can be imposed in a generic grid node and over border nodes of the area. Domain node coordinates are automatically generated by the code while a constant grid size, in each direction, was used. The code COVAR [45] was modified and enclosed into our code, to generate the statistical distribution of fracture apertures. COVAR uses the Choleski [23,30] matrix decomposition to produce the covariance matrix A A s2 exp ⫺ajlj;
14
where s 2 is the variance of Log apertures Y, a is the autocorrelation parameter and l is the separation ‘‘lag’’. The exponential form of the Eq. (14) suggests that we can define the correlation length equal to 2/a . The generation procedure of the apertures, Y log b, referred in each grid node, has been carried out using Y L1 ⫹ n;
15
where n is the Y mean value, e is a vector normally distributed with zero mean and standard deviation equal to 1 and L is the covariance decomposition matrix [11]. The standard deviation (0.3) and correlation length (1 m), for statistical generation of the fracture apertures, were derived by means of the results of field tracer tests carried out under natural (not pumped) fluid flow conditions [6,15]. The mean value may be estimated from the results of hydraulic and tracer tests, as proposed by Shapiro and Nicholas [37] or Ghelar [18]. In fact, using the residence time of the tracer, tr, the tracer aperture [42], btr , can be evaluated by 6m r 1=2 r2 ⫺ r02 log ; 16 btr gDftr r0 where r is the distance between the injection well and the central one, r0 is the pumping well radius and Df is the difference of piezometric heads. The mass balance aperture bb is given by bb
Qw × tr : p r 2 ⫺ r02
17
332
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Table 1 Hydraulic apertures (mm) in the East (E), West (W), North (N) and South (S) direction Qw (l/s) Df (m) bt 4:0 m bt 3:0 m Df (m) bt 4:0 m bt 3:0 m Df (m) bt 4:0 m bt 3:0 m
10.6
30.8
68.1
b b b b b b
E
W
N
S
0.058 172 200 0.168 175 204 0.478 154 187
0.038 213 246 0.128 201 234 0.398 169 205
0.028 248 287 0.158 181 211 0.538 146 176
0.021 267 331 0.121 206 241 0.386 172 208
The hydraulic aperture may be found using well drawdowns and solving the radial flow problem during pumping test [26], 6Qw m r 1=2 log : 18 b pDfbt g r0 Using a total aperture bt equal to bb , the Eq. (18) returns to the Eq. (16). By means of the results of hydraulic tests [14] some hydraulic apertures were estimated (Table 1) using two hypothetical values of the total aperture bt. These apertures may be very distant from the average value of fracture apertures. The effective aperture may be several times greater than tracer apertures when s of the fracture apertures is 0.3–0.4 and Peclet number is low [41]. For instance, using the Ghelar [18] equation, the hydraulic aperture of 0.270 mm should be increased to beff
b e
⫺1n 10s2 =2
1:5 × 270mm 艑 400mm
19
to obtain the mean value for the statistical fracture generation.
5. The code application The code was applied to simulate two-dimensional solute transport during pumping in a square area of a fractured aquifer. Because the length of the square border was 20 m, Dx, and consequently Dy, were selected equal to 1 m (441 nodes). Constant piezometric heads along all the border nodes and in the central one were used as model boundary conditions (Fig. 4) to reproduce the stationary cone of depression during pumping. Other effects, as for example vertical leakage, should be evaluated previously [25] and the code can be applied only when stationary conditions are achieved. The code solves the flow problems for each pumping discharge. System solutions, i.e. the piezometric heads, were found (Fig. 5). The tolerance in the SOR method was selected equal to 10 ⫺6 m. In the computational code each particle can be introduced in every grid node. Some characteristic particle situations were used to locate outlet particle conditions (as particles out of the square area or too slow). Several particles were left in the border nodes, to simulate the solute injection, and collected in the central one (pumping well, C). As known [44], the number of particles should be large enough so that the characteristics of the breakthrough curves do not change. The number of particles recovered in the outlet location was 1171, 1702 and 1890 during simulated pumping of respectively 10.6, 30.8 and 68.1 l/s. The number of particles left in the inlet locations was a few larger (2000). No tracer sorption has been supposed during computer simulations. The breakthrough histograms (BTH’s) (Fig. 7–8) were obtained by plotting the relative residence time frequency of the all particles collected in the outlet location (pumping well). An estimate of the BTH’s representativity can be realised placing on the top of BTH’s the analytical solutions obtained using fixed dispersion coefficients. In fact
Fig. 4. Boundary conditions.
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
333
Fig. 5. Piezometric water tables (in meters on sea water level) under steady state conditions before (a) and during pumping of 10.6 l/s (b), 30.8 l/s (c) and 68.1 l/s (d).
the same figure also contains the breakthrough curves (BTC’s) given by analytical solutions where the hydrodynamic dispersion coefficients were estimated by ‘‘channelling’’ approach [31]. Also if these BTC’s do not fit exactly the particle’s histograms (see Figs. 7–8), because the actual flow is not uniform and one-dimensional, as the analytical solution required, a sufficient approximation may be noted. This result is consistent with those achieved by some authors [4] for correspondent parameters, i.e. variance and integral scale, and dimensionless time ⬍ 2. Also the constancy of the variance with the mean distance covered from the particles was verified. The path of the each particle (Fig. 6) was recorded as well as the total residence time. By using the modal time of particles, tr (Table 2), the tracer aperture and the mass balance aperture were estimated by means of Eqs. (16) and (18). The effects of non-laminar flow on the solute transport,
were evaluated through a pumping simulation of Qw 400 l=s. The well drawdowns for this discharge were evaluated by the best-fit functions F w of the piezometric heads in each well
Fw C1 ⫹ C2 ·log Qw × Qw :
20
By using the C1 and C2 coefficients given by Table 3, determination coefficients of Eq. (20) ranging from 0.990– 0.997. The computational code can evaluate the friction factor by Eq. (8) and by imposing boundary conditions derived by Eq. (20), the piezometric heads under simulated pumping of 400 l/s were found. The estimated Reynolds numbers, in every grid channels showed non-laminar conditions also for Qw 30:8 l=s. The BTH’s obtained by code applications show different behaviours (Fig. 8) if respectively, laminar and non-laminar equation, is considered. The computed discharge in the central node of the grid
334
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Fig. 6. Pathways of the particles during pumping of 30.8 l/s.
model was 9.5 l/s, 30 l/s and 66 l/s when Q w was respectively 10.6 l/s, 30.8 l/s and 68.1 l/s. These values were obtained by Qw
p r2 ⫺ r20 ·bt tr
21
using a total aperture bt 3:0 m and the modal residence times of particles (see Table 2). The tracer residence time given by the model can be used to calibrate the mean aperture of the equivalent single fracture approach. The discrepancy between the calculated and observed discharges may be varied using different values of the mean aperture, for each discharge. The discharge errors evaluated, respectively ⫺10%, ⫺2% and ⫺3%, mean that the mean hydraulic aperture in the square area should be considered inferior when the pumping flow rate increases, as well as the trend of calculated hydraulic aperture shows (see Table. 1). An increase of the fluid-flow resistance in the square area could cause such a trend of the apertures, when the velocity, i.e. the Reynolds number, increases. Indeed the fluid flow resistance is related inversely to the hydraulic aperture (see Eq. (3) or Eq. (10)) and the latter decreases when the flow resistance increases.
6. Conclusion A computational code was applied to simulate solute transport under non-uniform and steady state flow conditions. The solution of Navier–Stokes and Darcy–Weisbach’s equation were used to solve the flow problem under variable flow’s regime. Some data of pumping tests were used to evaluate the mean aperture of the single fracture. This value was used for statistical generation of the apertures. The average aperture, estimated under different pump flow rates, seem to decrease when the discharge increases. An increasing of the fracture flow resistance may cause such trend, when the velocity, i.e. the Reynolds number, increases. Using the computational code results, the mass balance aperture and the tracer aperture were estimated. These values are consistent with hydraulic apertures, and showed similar trend. The simulated pumping results showed that under non-laminar conditions more accurate interpretations are required. The code may be useful to study pollutant transport in the capture zones of the well, to locate the real pollutant sources and draw control area of the wells for drinking uses. Nevertheless some limitations of the code application in vast areas can be found because the grid size cannot be very large.
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Fig. 7. Breakthrough histograms and curves during pumping (10.6 l/s (a), 30.8 l/s (b) and 68.1 l/s (c)).
335
336
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Fig. 8. Breakthrough histograms and curves during the pumping simulation of 400 l/s.
Table 2 Average aperture values Qw (l/s)
Df (m)
tr (d)
bb (m)
btr (mm)
10.6 30.8 68.1
0.03625 0.14375 0.45000
1.150 0.390 0.164
3.35 3.30 3.07
256 221 192
Table 3 Simulated drawdown for Q w 400 1=s in the Central (C), East (E), West (W), North (N) and South (S) well
C1(m) C2 (mm s/1) f (m)
C
E
W
N
S
1.404 ⫺ 0.767 ⫺ 0.434
1.421 ⫺ 0.438 ⫹ 0.370
1.394 ⫺ 0.476 ⫹ 0.254
1.377 ⫺ 0.357 ⫹ 0.522
1.379 ⫺ 0.469 ⫹ 0.256
Acknowledgements The author appreciated very much suggestions and comments given by C. Datei and A. Rinaldo of the Padova University and by M. Benedini of CNR-IRSA.
References [1] Abelin H, Birgerson L, Wide´n H, Agren T, Moreno L, e Neretnieks I. Channelling experiments in crystalline fractured rocks. Journal of Contaminant Hydrology 1994;15:129–158. [2] Abramovitz M, Stegun IA. In: Handbook of mathematical functions. New York: Dover, 1970, pp. 949–953.
[3] Bear J. Dynamics of fluids in porous media. New York: Elsevier, 1972, pp. 119–360. [4] Bellin A, Salandin P, Rinaldo A. Simulation of dispersion in heterogeneousporous formations: statistics first-order theories convergence of computations. Water Resources Research 1992;28:2211–2227. [5] Benedini M, Giuliano G, Troisi S. Experimental study of water flow in ractured media by means of a hydraulic model. In: Ganoulis J, Ministry of Civ. and of the Sciences, editors. Proc. int. symposium on scale effects in porous media. Thessaloniki GR, 1978, pp. 2.41– 2.55. [6] Benedini M, Di Fazio A, Masciopinto, C, Vurro, M. Interpretazione delle curve di restituzione di prove in situ per lo studio della dispersione in un sistema fessurato. Proc. XXIII conv. di idraulica e costruzioni idrauliche, Florence (I), August, 31-September. 4(1), 1992, pp. B15–B26. [7] Blumberg PN, Curl RL. Experimental and theoretical studies of dissolution roughness. J. Fluid Mech. 1974;65:735–751. [8] Brown SR. Fluid flow trough rock joint: the effect of surface roughness. Journal of Geophysical Research 1987;92:1337–1347. [9] Brusseau Mark L, Rao PSC. Sorption nonideaility during organic contaminant transport in porous media. Critical Reviews in Environmental Control 1989;19(I):33–99. [10] Bursik M. Theory of the sedimentation of suspended particles from fluvial plumes. Sedimentology 1995;42(6):831–838. [11] Clifton PM, Neuman SP. Effects of Kriging and inverse modelling on conditional simulation of the arva valley aquifer in Southern Arizona. Water Resources Research 1982;18(4):1215–1234. [12] Dagan G. Time-dependent macrodispersion for solute transport in anisotropic heterogeneous aquifers. Water Resources Research 1988;24:1491–1500. [13] Dagan G. Flow and transport in porous formations. New York: Springer-Verlag, 1989, pp. 157–350. [14] Di Fazio A, Maggiore M, Masciopinto C, Troisi S, Vurro M. Caratterizzazione idrogeologica di una stazione di misura in ambiente carsico. Acque Sotterranee 1992;2:17–26. [15] Di Fazio A, Masciopinto C, Troisi S, Vurro M, Una interpretazione bidimensionale di prove di tracciamento in sistemi fratturati. Proc.
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
[16]
[17] [18] [19]
[20]
[21] [22] [23]
[24] [25]
[26]
[27]
[28] [29] [30] [31]
XXIV conv. di idraulica e costruzioni idrauliche. Florence (I), September 20-22. vol. 1, 1993, pp. T1-167–T1-178. Di Fazio A, Masciopinto C, Troisi S. Tracer tests under non-linear laminar flow’s regime in a karstified aquifer, Proceedings of first international conference on the the impact of industry on groundwater resources, May 22–24, Cernobbio, Como, Italy, 1996, pp. 399–410. Gale SJ. The hydraulics of conduit flow in carbonate aquifer. J. Hydrol. 1984;70:309–327. Gelhar Lynn W. Stochastic subsurface hydrology, Englewood Cliffs, NJ: Prentice-Hall, 1993. Goode DJ. Particle velocity interpolation in block-centered finite difference grounwater flow models. Water Resouces Research 1990;26(5):925–940. Goode DJ, Shapiro AM. Comment on flow and tracer transport in a single fracture: a stochastic model and its relations to some field observations. Water Resource Research 1991;27(1):129–131. Gutfraind R, Hansen A. Study of fracture permeability using lattice gas automata. Transport in Porous Media 1995;18:131–149. Howard AD, Groves CG. Early development of karst systems 2. Turbulent flow, Water Resources Research, 1995;31:19–26. Larcombe MHE. In: Reid JK, editor. A list processing approach to the solution of large sparse sets of matrix equations and the factorisation of the overall matrix. Large sparse sets of linear equations. London: Academic Press, 1971, pp. 25–40. Marchi E, Rubatta A. Meccanica dei fluidi. UTET, Turin 1981;I:264– 448. Masciopinto C, Passarella G, Vurro M, Castellano L. Numerical simulations for the evaluation of the free surface history in porous media: comparison between two different approaches. Advances in Engineering Software 1994;21:149–157. Masciopinto C. L’utilizzo del particle tracking nello studio della propagazione degli inquinanti in sistemi carsici. IDROTECNICA 1995;4:245–253. Moreno L, Tsang YW, Tsang CF, Hale V, Neretnieks I. Flow and tracer transport in a single fracture: a stochastic model and its relation to some field observations. Water Resources Research 1988;24:2033–2048. Moreno L, Tsang YW, Tsang CF. Reply. Water Resources Research 1991;27(1):133–134. Moreno L, Neretnieks I. Fluid flow and solute transport in a network of channels. Journal of Contaminant Hydrology 1993;14:163–192. Nash JC. Compact numerical method for computers: linear algebra and function minimization. New York: Wiley, 1979. Neretnieks I. A note on fracture flow dispersion mechanisms in the ground. Water Resources Research 1983;19:364–370.
337
[32] Rinaldo A, Bellin A. Il ruolo delle eterogeneita` nell’analisi dei problemi di flusso e trasporto nei mezzi porosi. Proc. l’acqua nel sottosuolo: utilizzazione e salvaguardia. aspetti geologici, idraulici, geotecnici e ambientali. Taormina, May 11–13, vol. 4, 1995, pp. 37– 46. [33] Robinson PC. Connectivity of fracture systems-a percolation theory approach. J. Phys. A: Math. Gen. 1983;16:605–614. [34] Rose HE. An investigation into the laws of flow of fluids through beds of granular material. Proc. Inst. Mech. Eng., vol. 153, 1945, pp. 141– 148. [35] Salandin P, Rinaldo A, Dagan G. A note on transport in stratified formations by flow tilted with respect to the bedding. Water Resources Research 1991;27:3009–3017. [36] Schwartz F, Smith L, Crowe A. A stochastic analysis of macroscopic dispersion in fractured media. Water Resources Research 1983;19:1253–1265. [37] Shapiro AM, Nicholas JR. Assessing the validity of the channel model of fracture aperture under field conditions. Water Resources Research 1989;25:817–828. [38] Snow DT. The frequency of aperture of fracture in rock. Int. Journal of Rock Mechanics and Mineral Science 1970;7:23–40. [39] Tsang YW, Witherspoon PA. Hydromechanical behaviour of a deformable rock fracture subject to normal stress. Journal of Geophysical Research 1981;86:9287–9298. [40] Tsang YW, Tsang CF. Channel model of flow through fractured media. Water Resources Research 1987;23:467–479. [41] Tsang YW, Tsang CF, Neretnieks I, Moreno L. Flow and tracer transport in fractured media:a variable aperture channel model and its properties. Water Resources Research 1988;24:2049–2060. [42] Tsang YW. Usage of equivalent apertures for rock fractures as derived from hydraulic and tracer tests. Water Resources Research 1992;28:1451–1455. [43] Vandergraaf TT, Park CK, Drew DJ. Migration of conservative and poorly sorbing tracers in granite fractures. Proc. Fifth Int. Conf. On High Level Waste Management, Las Vegas, NV, USA: American Nuclear Soc., 1994, pp. 2251–2256. [44] Wille Nordqvist A, Tsang YW, Tsang CF, Bjo¨rn D, Andersson J. A variable aperture fracture network model for flow and transport in fractured rocks. Water Resources Research 1992;28: 1703–1713. [45] Williams SA, El Kadi AI. COVAR- A computer program for generating two-dimensional fields of autocorrelated parameters by matrix decomposition. IGWMC -International groundwater modelling center, c/o tno Deft, Netherlands: Institute of Applied Geoscience, 1986.
Advances in Engineering Software 30 (1999) 327–337
Particles’ transport in a single fracture under variable flow regimes C. Masciopinto* National Research Council, Water Research Institute, Via F. De Blasio 5, 70123 Bari, Italy Received 1 June 1996; received in revised form 1 October 1998; accepted 1 October 1998
Abstract A mathematical model to account for non-laminar terms, during simulated transports in fractured rocks has been developed. The relevant conceptual model consists of a two dimensional grid block equivalent to a variable fracture plane. The block aperture variability is expressed by means of stochastic procedures, based upon some characteristics of fractured media. One of them, e.g. the equivalent hydraulic aperture, may be estimated by means of the results of field tests. The solution of Navier–Stokes and Darcy–Weisbach’s equation for steady state flow between two parallel plates with constant aperture allows the systems of equations under respetively, laminar and non-laminar conditions, to be wirtten. The flow model is a standard finite-difference method with a modified version of Darcy’s equation to account for non-laminar terms that are important at high water velocity. To simulate solute transports, the particle tracking technique was applied. 䉷 1999 Elsevier Science Ltd. All rights reserved. Keywords: Fractured media; Radial flow; Non-laminar regime; Numerical model; Fracture apertures
1. Introduction Particle tracking technique can be very useful for the identification of the control area in well field. In this article such a technique is used to study the tracer transport in a square area of a fractured aquifer. The study of pollutant transport in an aquifer cannot neglect the physical characteristics of soils, being either fractured or granular. The usual diffusive model, which describes flow and transport of pollutants in terms of an ‘‘equivalent porous media’’, may be inconsistent in fractured systems. In fact, in systems with few fractures or in a heterogeneous porous media, the basic assumption of the continuum approach, i.e. the Representative Elementary Volume (REV) to which all physical properties of the medium are referred, is often missing. Recently the efficacy of the stochastic approach to study heterogeneous porous formations has been showed [4,35]. Medium heterogeneity can be described in a schematic representation, also three-dimensional, by means of significant parameters [12]. In rock systems with few fractures some authors [8,21] proposed fractal theory, to reproduce a single fracture and, to study the fluid flow; the particle tracking can be applied, as clearly explained by Schwartz et al. [36], to study pollutants transport in 2- or 3-D scale [1,29,42,44]. This method * Corresponding author. Fax: 0039 80 5313365. E-mail address: [email protected] (C. Masciopinto)
was applied by Schwartz et al. [36] to study solute transport in fracture sets having constant apertures and different orientations. For Neretnieks [31] and Tsang and Tsang [40], a single fracture can be considered equivalent to a set of independent channels, with variable cross-section. The apertures can be derived by log-normal distributions [38] of values and can be generated by means of an exponential function of the spatial covariance of the apertures. Fracture apertures can be obtained when the values of mean, standard deviation and correlation length are known. Moreno et al. [27] applied this technique to study the dispersion through a single fracture with variable aperture under uniform and laminar one-dimensional flow conditions. A computational code has been developed here, to model non-uniform steady state flow under variable flow regimes in a variable-aperture fracture plane (Fig. 1). These groundwater fluid flow conditions are frequently verified inside the capture zones of pumping wells. The flow model is based on a standard finite difference method with a modified version of Darcy’s equation, to account for non-linear terms that are important at high velocity values. An application of this model allows the evaluation of the non-laminar effects on the pollutant transport in fractures to be done. Clogging effects and related phenomena have not been considered in this article. The fact that the fracture aperture is as low as 1–20 mm, in the narrowest parts of many fracture channels, means that most of the soil particles that can be moved, belong to the clay fraction. Clay minerals tend to
0965-9978/99/$ - see front matter 䉷 1999 Elsevier Science Ltd. All rights reserved. PII: S0965-997 8(98)00092-1
328
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Fig. 1. Fracture generation parameters for simplified models.
be stuck by interacting with fracture minerals. Thus, the drug forces need to overcome bonds between transported soil particles and fracture minerals. To avoid clogging effects caused by transported soil particles a flow rate exceeding a minimum value (about 8.6 m/d [10]) is required. These effects will be considered in the next article.
2. Fluid flow motion Fluid flow can be studied applying the solution of the Navier–Stokes’ equation for average velocity in a single fracture of constant width, b, bounded by two parallel impervious plates [39]
g 2 b 7f; V ⫺ 12m
1
where m and g are respectively the fluid dynamic viscosity and specific weight (g /m is about 10 7 m ⫺1s ⫺1 for water at
20⬚C), while f is the piezometric head in the extremities of the channel with constant length Dx and constant crosssection Dy × b. Using a finite difference method with a block centred approach, the grid nodes (Fig. 2) may be defined using a grid discretization step size capable to replace flow streamlines with subsequent segments in x and y direction. Moreover, if i and j are two consecutive grid nodes, oriented channels in x or y direction can be defined. As the flow is parallel to the grid, the effects of grid orientation on the flow velocity [19] are not significant [28]. These channels can be considered with constant length Dx , or Dy, and constant cross-section Dy × b, or Dx × b. Discharge values can be obtained from the equation Qij
g b3 Df Dy Dx m 12
2
Dx and Dy were selected short enough so that one-dimensional flow over each channel occurred and the relation
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
329
Fig. 2. Grid nodes and channels.
between interconnected piezometric heads and discharge is given by
fi ⫺ fj
Dfi ⫹ Dfj Qij 2 1=12 g=mb3i Dy 2=Dx Qij 1=12 g=mb3j Dy 2=Dx !# 1 1 ⫹ 3 Qij Rij ; b3i bj
Qij 6
m Dx g Dy
X
Qij
J
⫹
"
grid, without sink and source, we can write
3
where Rij is the flow resistance of the channel between nodes i and j under laminar regime. Eq. (3) proposed a linear gradient interpolation between adjacent blocks with constant apertures bi and bj. In fact, using bi and bj into Eq. (2), two piezometric head differences can be found, Df i and Df j, respectively. Thus, because of nodal step changes in thickness, the linear interpolation of the gradient gives a channel conductance, 1/Rij, equal to the harmonic mean of the nodal values. To preserve the net discharge across each node of the
X fi ⫺ fj 0; Rij J
4
where the converging channels, J, are three, instead of four, along the borders of the studied area. By using Eq. (4) for all grid nodes a system of linear equations can be obtained; the system solution can be found after imposing suitable boundary conditions. To improve this method, the flow resistances were corrected to consider inertial forces for high Reynolds numbers (⬎ 10). Under these field conditions, the Fanning factor [3] f can be used to evaluate the frictional losses inside each channel f
2g Df ; V 2 Dx
5
where g is the acceleration of gravity. The Fanning friction factor in laminar condition, for fluid flow motion between two parallel impervious surfaces, can be derived combining
330
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Fig. 3. Relationships between the Fanning friction factor and Reynolds number.
Eqs. (1) and (5). If the hydraulic radius Ri is Ri
Dy × b b 艑 ; 2Dy ⫹ 2b 2
6
valid for b ⬍⬍ Dy: Fanning values can be obtained as f 24=Re, where Re
4rVRi b 2rV m m
7
is the Reynolds number for channel with constant rectangular cross-section [24]. Also for laminar flow through porous (granular) media, the relationship between Fanning friction factor and Reynolds number gives a straight line portion, in a log–log plane, expressed by f 1000=Re. Thus under laminar conditions the f f Refunctions, are two parallel straight lines in a log–log plane for both porous media and impervious surfaces. Under non-laminar conditions, in agreement with Gale [17], measurements of the friction factors are rare, despite the importance of such values in hydrological studies. Results obtained by Blumberg and Curl [7] and Benedini et al. [5], are not suitable in case of limestone that does not have very large conduits (caves) or regular channels (e.g., bricks). Assuming that the shape of relations f f Re holds also for non-laminar regimes, the Fanning friction factor for flow through a single fracture can be obtained by a translation of the granular porous media data, using the constant 24/1000 (Fig. 3). In this way the Darcy–Weisbach friction factor, l 4·f ; can be evaluated by means of the best-fit of porous media data [34], shifted by value 96/1000 105:12 992:31 300:70 96 ⫺ : 8 l 15:53 ⫹ p ⫹ Re 1000 eRe Re This equation gives a best fit also for laminar (Re ⬍ 2) and
turbulent [22] (Re ⬎ 10000) regimes. To verify this behaviour, a forced tracer test with high discharges was carried out and some points under non-laminar regime were found [16] (see Fig. 3). Using the Darcy–Weisbach equation, the corrected discharge can be obtained by s 4bg Df Qij b × Dy 9 × l Dx after that, the solution of the system (4) was found and the Re numbers calculated by means of the velocity values given by Eq. (1). The relation between interconnected piezometric heads and discharge, under generic flow regime, is now given by ! Dfi ⫹ Dfj Q2ij li;j 1 1 × fi ⫺ fj ⫹ 3 2 8gDy Dy=Dx b3i bj " !# l i;j Dx 1 1 × ⫹ 3 Q ij 2 D y b3i 8g D y bj Q ij
2
R 0ij ; 10
Rij0
is the corrected flow resistance of the channel where between nodes i and j. Without sinks and sources, in each grid node, we can write s X fi ⫺ fj X 0 11 Qij R 0ij that represents a system of non-linear equations, the solution of which can be found after imposing suitable boundary conditions. The developed computational code employs the successive over relaxation method (SOR) to solve
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
systems (4) and (11), after linearization. When laminar conditions occur, systems (4) and (11) give the same solution.
3. Tracer transport The phenomenon of solute transport is investigated by tracking the particles advected through an ideal single fracture with variable apertures and with a behaviour which is equivalent to that of a fractured media (see Fig. 1). The particles can be introduced into a generic node of the domain. After that, the constant discharges within each channel were evaluated. However, during pumping simulations, the particles have been left in border nodes and collected in the central node of the domain. Particles coming to an intersection (i.e., the generic block-centred grid node) are distributed in the outlet branches with a probability proportional to the flow rate [2,33]. This method is consistent with an instantaneous and complete solute mixing in each node and the majority of particles take the flow path with higher discharge. If dM nC0 x0 dx) is the polluter mass with concentration C0 at t t0 at the inlet of the channel, its Lagrangian coordinates can be expressed as a function of the probability density of the flow path particles [13,32]. Here, the numerical solution was used. Following each particle along its pathway as the Lagrangian point of view requires, it migrates inside channels, oriented in the x or y direction (see Fig. 2), consecutively. Within each channel uniform and steady-state flow occurs and for preserving the discharge, as the channel cross section is constant, also the velocity must be constant (incompressible fluid and non-deformable rock). To observe these hypotheses the size of the grid block cannot be very large. As the mass flow at outlet location is proportional to the number of particles that reach that location, the concentration is derived from the number of particles per unit time that arrive at a collection point divided by the flowrate at the same location [29]. The mean residence time of the particle inside each grid element of constant length Dx and constant cross-section Dy × bij (i.e., channels in x direction) can be determined by using Dx ti ; Vij
12
where Vij is given by arithmetic mean of the apertures (see Fig. 2). As Dx Dy the same equation can be applied for channels oriented in the y direction. This equation was preferred to that proposed by Moreno et al. [27] to avoid lateral effects on the volumetric flux from node i to node j [20]. The total particle residence time is given by the incremental addition of these times over the entire particle path from the inlet to the outlet position. Following Moreno et al. [27], under equilibrium conditions, a sorbing tracer in each channel can be considered
331
using the retardation factor, Ra, for each node Ra 1 ⫹ 2
Ka ; bi
13
where Ka is the surface distribution coefficient which can be obtained from static sorption experiments (i.e., batch tests) [9,43]. In fractured systems with large macropores, the tracer diffusion into fracture walls cannot be neglected and kinetic equations should be adopted to consider the effects of diffusion. 4. The fracture generation The computational code has been written to solve fluid flow and pollutant transport in the plane of the statistically generated fracture. To solve the fluid flow problem, steady state boundary conditions can be imposed in a generic grid node and over border nodes of the area. Domain node coordinates are automatically generated by the code while a constant grid size, in each direction, was used. The code COVAR [45] was modified and enclosed into our code, to generate the statistical distribution of fracture apertures. COVAR uses the Choleski [23,30] matrix decomposition to produce the covariance matrix A A s2 exp ⫺ajlj;
14
where s 2 is the variance of Log apertures Y, a is the autocorrelation parameter and l is the separation ‘‘lag’’. The exponential form of the Eq. (14) suggests that we can define the correlation length equal to 2/a . The generation procedure of the apertures, Y log b, referred in each grid node, has been carried out using Y L1 ⫹ n;
15
where n is the Y mean value, e is a vector normally distributed with zero mean and standard deviation equal to 1 and L is the covariance decomposition matrix [11]. The standard deviation (0.3) and correlation length (1 m), for statistical generation of the fracture apertures, were derived by means of the results of field tracer tests carried out under natural (not pumped) fluid flow conditions [6,15]. The mean value may be estimated from the results of hydraulic and tracer tests, as proposed by Shapiro and Nicholas [37] or Ghelar [18]. In fact, using the residence time of the tracer, tr, the tracer aperture [42], btr , can be evaluated by 6m r 1=2 r2 ⫺ r02 log ; 16 btr gDftr r0 where r is the distance between the injection well and the central one, r0 is the pumping well radius and Df is the difference of piezometric heads. The mass balance aperture bb is given by bb
Qw × tr : p r 2 ⫺ r02
17
332
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Table 1 Hydraulic apertures (mm) in the East (E), West (W), North (N) and South (S) direction Qw (l/s) Df (m) bt 4:0 m bt 3:0 m Df (m) bt 4:0 m bt 3:0 m Df (m) bt 4:0 m bt 3:0 m
10.6
30.8
68.1
b b b b b b
E
W
N
S
0.058 172 200 0.168 175 204 0.478 154 187
0.038 213 246 0.128 201 234 0.398 169 205
0.028 248 287 0.158 181 211 0.538 146 176
0.021 267 331 0.121 206 241 0.386 172 208
The hydraulic aperture may be found using well drawdowns and solving the radial flow problem during pumping test [26], 6Qw m r 1=2 log : 18 b pDfbt g r0 Using a total aperture bt equal to bb , the Eq. (18) returns to the Eq. (16). By means of the results of hydraulic tests [14] some hydraulic apertures were estimated (Table 1) using two hypothetical values of the total aperture bt. These apertures may be very distant from the average value of fracture apertures. The effective aperture may be several times greater than tracer apertures when s of the fracture apertures is 0.3–0.4 and Peclet number is low [41]. For instance, using the Ghelar [18] equation, the hydraulic aperture of 0.270 mm should be increased to beff
b e
⫺1n 10s2 =2
1:5 × 270mm 艑 400mm
19
to obtain the mean value for the statistical fracture generation.
5. The code application The code was applied to simulate two-dimensional solute transport during pumping in a square area of a fractured aquifer. Because the length of the square border was 20 m, Dx, and consequently Dy, were selected equal to 1 m (441 nodes). Constant piezometric heads along all the border nodes and in the central one were used as model boundary conditions (Fig. 4) to reproduce the stationary cone of depression during pumping. Other effects, as for example vertical leakage, should be evaluated previously [25] and the code can be applied only when stationary conditions are achieved. The code solves the flow problems for each pumping discharge. System solutions, i.e. the piezometric heads, were found (Fig. 5). The tolerance in the SOR method was selected equal to 10 ⫺6 m. In the computational code each particle can be introduced in every grid node. Some characteristic particle situations were used to locate outlet particle conditions (as particles out of the square area or too slow). Several particles were left in the border nodes, to simulate the solute injection, and collected in the central one (pumping well, C). As known [44], the number of particles should be large enough so that the characteristics of the breakthrough curves do not change. The number of particles recovered in the outlet location was 1171, 1702 and 1890 during simulated pumping of respectively 10.6, 30.8 and 68.1 l/s. The number of particles left in the inlet locations was a few larger (2000). No tracer sorption has been supposed during computer simulations. The breakthrough histograms (BTH’s) (Fig. 7–8) were obtained by plotting the relative residence time frequency of the all particles collected in the outlet location (pumping well). An estimate of the BTH’s representativity can be realised placing on the top of BTH’s the analytical solutions obtained using fixed dispersion coefficients. In fact
Fig. 4. Boundary conditions.
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
333
Fig. 5. Piezometric water tables (in meters on sea water level) under steady state conditions before (a) and during pumping of 10.6 l/s (b), 30.8 l/s (c) and 68.1 l/s (d).
the same figure also contains the breakthrough curves (BTC’s) given by analytical solutions where the hydrodynamic dispersion coefficients were estimated by ‘‘channelling’’ approach [31]. Also if these BTC’s do not fit exactly the particle’s histograms (see Figs. 7–8), because the actual flow is not uniform and one-dimensional, as the analytical solution required, a sufficient approximation may be noted. This result is consistent with those achieved by some authors [4] for correspondent parameters, i.e. variance and integral scale, and dimensionless time ⬍ 2. Also the constancy of the variance with the mean distance covered from the particles was verified. The path of the each particle (Fig. 6) was recorded as well as the total residence time. By using the modal time of particles, tr (Table 2), the tracer aperture and the mass balance aperture were estimated by means of Eqs. (16) and (18). The effects of non-laminar flow on the solute transport,
were evaluated through a pumping simulation of Qw 400 l=s. The well drawdowns for this discharge were evaluated by the best-fit functions F w of the piezometric heads in each well
Fw C1 ⫹ C2 ·log Qw × Qw :
20
By using the C1 and C2 coefficients given by Table 3, determination coefficients of Eq. (20) ranging from 0.990– 0.997. The computational code can evaluate the friction factor by Eq. (8) and by imposing boundary conditions derived by Eq. (20), the piezometric heads under simulated pumping of 400 l/s were found. The estimated Reynolds numbers, in every grid channels showed non-laminar conditions also for Qw 30:8 l=s. The BTH’s obtained by code applications show different behaviours (Fig. 8) if respectively, laminar and non-laminar equation, is considered. The computed discharge in the central node of the grid
334
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Fig. 6. Pathways of the particles during pumping of 30.8 l/s.
model was 9.5 l/s, 30 l/s and 66 l/s when Q w was respectively 10.6 l/s, 30.8 l/s and 68.1 l/s. These values were obtained by Qw
p r2 ⫺ r20 ·bt tr
21
using a total aperture bt 3:0 m and the modal residence times of particles (see Table 2). The tracer residence time given by the model can be used to calibrate the mean aperture of the equivalent single fracture approach. The discrepancy between the calculated and observed discharges may be varied using different values of the mean aperture, for each discharge. The discharge errors evaluated, respectively ⫺10%, ⫺2% and ⫺3%, mean that the mean hydraulic aperture in the square area should be considered inferior when the pumping flow rate increases, as well as the trend of calculated hydraulic aperture shows (see Table. 1). An increase of the fluid-flow resistance in the square area could cause such a trend of the apertures, when the velocity, i.e. the Reynolds number, increases. Indeed the fluid flow resistance is related inversely to the hydraulic aperture (see Eq. (3) or Eq. (10)) and the latter decreases when the flow resistance increases.
6. Conclusion A computational code was applied to simulate solute transport under non-uniform and steady state flow conditions. The solution of Navier–Stokes and Darcy–Weisbach’s equation were used to solve the flow problem under variable flow’s regime. Some data of pumping tests were used to evaluate the mean aperture of the single fracture. This value was used for statistical generation of the apertures. The average aperture, estimated under different pump flow rates, seem to decrease when the discharge increases. An increasing of the fracture flow resistance may cause such trend, when the velocity, i.e. the Reynolds number, increases. Using the computational code results, the mass balance aperture and the tracer aperture were estimated. These values are consistent with hydraulic apertures, and showed similar trend. The simulated pumping results showed that under non-laminar conditions more accurate interpretations are required. The code may be useful to study pollutant transport in the capture zones of the well, to locate the real pollutant sources and draw control area of the wells for drinking uses. Nevertheless some limitations of the code application in vast areas can be found because the grid size cannot be very large.
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Fig. 7. Breakthrough histograms and curves during pumping (10.6 l/s (a), 30.8 l/s (b) and 68.1 l/s (c)).
335
336
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
Fig. 8. Breakthrough histograms and curves during the pumping simulation of 400 l/s.
Table 2 Average aperture values Qw (l/s)
Df (m)
tr (d)
bb (m)
btr (mm)
10.6 30.8 68.1
0.03625 0.14375 0.45000
1.150 0.390 0.164
3.35 3.30 3.07
256 221 192
Table 3 Simulated drawdown for Q w 400 1=s in the Central (C), East (E), West (W), North (N) and South (S) well
C1(m) C2 (mm s/1) f (m)
C
E
W
N
S
1.404 ⫺ 0.767 ⫺ 0.434
1.421 ⫺ 0.438 ⫹ 0.370
1.394 ⫺ 0.476 ⫹ 0.254
1.377 ⫺ 0.357 ⫹ 0.522
1.379 ⫺ 0.469 ⫹ 0.256
Acknowledgements The author appreciated very much suggestions and comments given by C. Datei and A. Rinaldo of the Padova University and by M. Benedini of CNR-IRSA.
References [1] Abelin H, Birgerson L, Wide´n H, Agren T, Moreno L, e Neretnieks I. Channelling experiments in crystalline fractured rocks. Journal of Contaminant Hydrology 1994;15:129–158. [2] Abramovitz M, Stegun IA. In: Handbook of mathematical functions. New York: Dover, 1970, pp. 949–953.
[3] Bear J. Dynamics of fluids in porous media. New York: Elsevier, 1972, pp. 119–360. [4] Bellin A, Salandin P, Rinaldo A. Simulation of dispersion in heterogeneousporous formations: statistics first-order theories convergence of computations. Water Resources Research 1992;28:2211–2227. [5] Benedini M, Giuliano G, Troisi S. Experimental study of water flow in ractured media by means of a hydraulic model. In: Ganoulis J, Ministry of Civ. and of the Sciences, editors. Proc. int. symposium on scale effects in porous media. Thessaloniki GR, 1978, pp. 2.41– 2.55. [6] Benedini M, Di Fazio A, Masciopinto, C, Vurro, M. Interpretazione delle curve di restituzione di prove in situ per lo studio della dispersione in un sistema fessurato. Proc. XXIII conv. di idraulica e costruzioni idrauliche, Florence (I), August, 31-September. 4(1), 1992, pp. B15–B26. [7] Blumberg PN, Curl RL. Experimental and theoretical studies of dissolution roughness. J. Fluid Mech. 1974;65:735–751. [8] Brown SR. Fluid flow trough rock joint: the effect of surface roughness. Journal of Geophysical Research 1987;92:1337–1347. [9] Brusseau Mark L, Rao PSC. Sorption nonideaility during organic contaminant transport in porous media. Critical Reviews in Environmental Control 1989;19(I):33–99. [10] Bursik M. Theory of the sedimentation of suspended particles from fluvial plumes. Sedimentology 1995;42(6):831–838. [11] Clifton PM, Neuman SP. Effects of Kriging and inverse modelling on conditional simulation of the arva valley aquifer in Southern Arizona. Water Resources Research 1982;18(4):1215–1234. [12] Dagan G. Time-dependent macrodispersion for solute transport in anisotropic heterogeneous aquifers. Water Resources Research 1988;24:1491–1500. [13] Dagan G. Flow and transport in porous formations. New York: Springer-Verlag, 1989, pp. 157–350. [14] Di Fazio A, Maggiore M, Masciopinto C, Troisi S, Vurro M. Caratterizzazione idrogeologica di una stazione di misura in ambiente carsico. Acque Sotterranee 1992;2:17–26. [15] Di Fazio A, Masciopinto C, Troisi S, Vurro M, Una interpretazione bidimensionale di prove di tracciamento in sistemi fratturati. Proc.
C. Masciopinto / Advances in Engineering Software 30 (1999) 327–337
[16]
[17] [18] [19]
[20]
[21] [22] [23]
[24] [25]
[26]
[27]
[28] [29] [30] [31]
XXIV conv. di idraulica e costruzioni idrauliche. Florence (I), September 20-22. vol. 1, 1993, pp. T1-167–T1-178. Di Fazio A, Masciopinto C, Troisi S. Tracer tests under non-linear laminar flow’s regime in a karstified aquifer, Proceedings of first international conference on the the impact of industry on groundwater resources, May 22–24, Cernobbio, Como, Italy, 1996, pp. 399–410. Gale SJ. The hydraulics of conduit flow in carbonate aquifer. J. Hydrol. 1984;70:309–327. Gelhar Lynn W. Stochastic subsurface hydrology, Englewood Cliffs, NJ: Prentice-Hall, 1993. Goode DJ. Particle velocity interpolation in block-centered finite difference grounwater flow models. Water Resouces Research 1990;26(5):925–940. Goode DJ, Shapiro AM. Comment on flow and tracer transport in a single fracture: a stochastic model and its relations to some field observations. Water Resource Research 1991;27(1):129–131. Gutfraind R, Hansen A. Study of fracture permeability using lattice gas automata. Transport in Porous Media 1995;18:131–149. Howard AD, Groves CG. Early development of karst systems 2. Turbulent flow, Water Resources Research, 1995;31:19–26. Larcombe MHE. In: Reid JK, editor. A list processing approach to the solution of large sparse sets of matrix equations and the factorisation of the overall matrix. Large sparse sets of linear equations. London: Academic Press, 1971, pp. 25–40. Marchi E, Rubatta A. Meccanica dei fluidi. UTET, Turin 1981;I:264– 448. Masciopinto C, Passarella G, Vurro M, Castellano L. Numerical simulations for the evaluation of the free surface history in porous media: comparison between two different approaches. Advances in Engineering Software 1994;21:149–157. Masciopinto C. L’utilizzo del particle tracking nello studio della propagazione degli inquinanti in sistemi carsici. IDROTECNICA 1995;4:245–253. Moreno L, Tsang YW, Tsang CF, Hale V, Neretnieks I. Flow and tracer transport in a single fracture: a stochastic model and its relation to some field observations. Water Resources Research 1988;24:2033–2048. Moreno L, Tsang YW, Tsang CF. Reply. Water Resources Research 1991;27(1):133–134. Moreno L, Neretnieks I. Fluid flow and solute transport in a network of channels. Journal of Contaminant Hydrology 1993;14:163–192. Nash JC. Compact numerical method for computers: linear algebra and function minimization. New York: Wiley, 1979. Neretnieks I. A note on fracture flow dispersion mechanisms in the ground. Water Resources Research 1983;19:364–370.
337
[32] Rinaldo A, Bellin A. Il ruolo delle eterogeneita` nell’analisi dei problemi di flusso e trasporto nei mezzi porosi. Proc. l’acqua nel sottosuolo: utilizzazione e salvaguardia. aspetti geologici, idraulici, geotecnici e ambientali. Taormina, May 11–13, vol. 4, 1995, pp. 37– 46. [33] Robinson PC. Connectivity of fracture systems-a percolation theory approach. J. Phys. A: Math. Gen. 1983;16:605–614. [34] Rose HE. An investigation into the laws of flow of fluids through beds of granular material. Proc. Inst. Mech. Eng., vol. 153, 1945, pp. 141– 148. [35] Salandin P, Rinaldo A, Dagan G. A note on transport in stratified formations by flow tilted with respect to the bedding. Water Resources Research 1991;27:3009–3017. [36] Schwartz F, Smith L, Crowe A. A stochastic analysis of macroscopic dispersion in fractured media. Water Resources Research 1983;19:1253–1265. [37] Shapiro AM, Nicholas JR. Assessing the validity of the channel model of fracture aperture under field conditions. Water Resources Research 1989;25:817–828. [38] Snow DT. The frequency of aperture of fracture in rock. Int. Journal of Rock Mechanics and Mineral Science 1970;7:23–40. [39] Tsang YW, Witherspoon PA. Hydromechanical behaviour of a deformable rock fracture subject to normal stress. Journal of Geophysical Research 1981;86:9287–9298. [40] Tsang YW, Tsang CF. Channel model of flow through fractured media. Water Resources Research 1987;23:467–479. [41] Tsang YW, Tsang CF, Neretnieks I, Moreno L. Flow and tracer transport in fractured media:a variable aperture channel model and its properties. Water Resources Research 1988;24:2049–2060. [42] Tsang YW. Usage of equivalent apertures for rock fractures as derived from hydraulic and tracer tests. Water Resources Research 1992;28:1451–1455. [43] Vandergraaf TT, Park CK, Drew DJ. Migration of conservative and poorly sorbing tracers in granite fractures. Proc. Fifth Int. Conf. On High Level Waste Management, Las Vegas, NV, USA: American Nuclear Soc., 1994, pp. 2251–2256. [44] Wille Nordqvist A, Tsang YW, Tsang CF, Bjo¨rn D, Andersson J. A variable aperture fracture network model for flow and transport in fractured rocks. Water Resources Research 1992;28: 1703–1713. [45] Williams SA, El Kadi AI. COVAR- A computer program for generating two-dimensional fields of autocorrelated parameters by matrix decomposition. IGWMC -International groundwater modelling center, c/o tno Deft, Netherlands: Institute of Applied Geoscience, 1986.