Permeameter Data Verify New Turbulence Process for MODFLOW

Methods Note/

Permeameter Data Verify New Turbulence Process for MODFLOW by Eve L. Kuniansky1, Keith J. Halford2, and W. Barclay Shoemaker3

Abstract A sample of Key Largo Limestone from southern Florida exhibited turbulent flow behavior along three orthogonal axes as reported in recently published permeameter experiments. The limestone sample was a cube measuring 0.2 m on edge. The published nonlinear relation between hydraulic gradient and discharge was simulated using the turbulent flow approximation applied in the Conduit Flow Process (CFP) for MODFLOW-2005 mode 2, CFPM2. The good agreement between the experimental data and the simulated results verifies the utility of the approach used to simulate the effects of turbulent flow on head distributions and flux in the CFPM2 module of MODFLOW-2005.

Introduction Bulk density, total porosity, effective porosity, constant head, and permeameter test data were published for 13 cubes of Key Largo Limestone from southern Florida (DiFrenna et al. 2007). Seventy percent of the permeameter tests conducted on these samples remained under laminar flow, even with the extreme gradients imposed during the permeameter experiments (tables 2 and 3, DiFrenna et al. 2007). Published data along three orthogonal axes of flow with cube number 6 exhibited a nonlinear relation between hydraulic gradient and discharge (DiFrenna et al. 2007, figure 6). This nonlinear relation is indicative of turbulent or non-Darcian flow. Cube number 6 measures 0.2 m on an edge with a bulk density of 1.38 g/cc, a total porosity of 0.49, an effective porosity of 0.34, and a representative pore diameter of 0.01 m (DiFrenna et al. 2007). These data provide a data set from which to evaluate the utility of a simple, turbulent flow approximation recently released for MODFLOW-2005 (Shoemaker et al. 2008). Similar to pipe flow, laminar and turbulent conditions in porous media are differentiated by the value of

1Corresponding author: U.S. Geological Survey, 3850 Holcomb Bridge Road, Suite 160, Norcross, GA 30092; (770) 409-7716. 2U.S. Geological Survey, 2730 N. Deer Run Road, Carson City, NV 89701. 3U.S. Geological Survey, 3110 SW 9th Avenue, Ft. Lauderdale, FL 33315. Received November 2007, accepted March 2008. Journal compilation ª 2008 National Ground Water Association. No claim to original US government works. doi: 10.1111/j.1745-6584.2008.00458.x

768

the Reynolds number (Darcy 1856; Schneebeli 1955). The Reynolds number (Re) is the ratio of inertial and viscous forces and for porous media is defined as: qvd Re ¼ ð1Þ l where q is density of water (997.5 kg/m3 at 23.5
C), v is specific velocity or mean velocity and is calculated from total discharge divided by cross-sectional area (m/s), d is mean pore diameter (m), and l is the dynamic or absolute viscosity of water (0.000924 kg/ms at 23.5
C). Flow in porous media typically is laminar under natural gradients because pore diameters are small (less than 0.01 m). As a result of the small hydraulic radius (cross-sectional area divided by wetted perimeter), surface tension and frictional forces at the pore walls slow velocity; thus viscous forces dominate the inertial forces and flow remains laminar. The critical Reynolds number (NRe) defines the threshold where laminar flow ceases and turbulent flow begins. The term turbulent flow includes non-Darcian laminar flow in this study because the relation between hydraulic gradient and specific velocity is nonlinear. Typically, NRe ranges between 1 and 60 in porous media (Bear 1979; Freeze and Cherry 1979; Darcy 1856; Schneebeli 1955). The NRe for samples of different aquifer or experimental material can be determined through laboratory experiments and is determined from the mean velocity at the point where the hydraulic gradient and specific velocity cease to be proportional (Darcy 1856; Schneebeli 1955). Functional relations between turbulent head loss and rock properties are not well defined, as is the case for pipe flow (Geankoplis 1993). Turbulent head losses in

Vol. 46, No. 5—GROUND WATER—September–October 2008 (pages 768–771)

pumped wells and flowmeters have been simulated previously in MODFLOW using pipe-flow friction factor relations (Halford 2000). The empirical Sacham formula (Geankoplis 1993) worked well for this application because turbulent flow predominantly occurs in pipes. These empirical formulas derived for pipe friction factors can behave poorly in porous media because the NRe are 2 to 3 orders of magnitude less than in pipes. Turbulent head loss in porous media can be well described with a simple power-relation function of the Re rather than complex pipe-flow equations.

Turbulent Flow in MODFLOW The Conduit Flow Process (CFP) is a recently released module for MODFLOW-2005 (Harbaugh 2005) that can simulate turbulent flow in karstic aquifers (Shoemaker et al. 2008). The CFP simulates turbulent flow as a discrete pipe network connected to the MODFLOW cells via head-dependent flux terms, mode 1 (CFPM1) or as continuous preferential flow layers where flow can transition between laminar and turbulent in the row and column direction, mode 2 (CFPM2). Additionally, both a pipe network and preferential flow layers can be simulated in the same model, mode 3 (CFPM3). Conduit flow layers in CFPM2 simulate preferential flow through interconnected pores and vugs more than 0.01 m in diameter, such as those that are present in the carbonate rock making up the Biscayne Aquifer in southern Florida. CFPM2 requires less data than CFPM1 or CFPM3 because discrete pipe networks are not simulated. However, CFPM1 is designed to simulate laminar and turbulent flow in large networks of submerged conduits. CFPM1 is an update of subroutines for MODFLOW-2005 from the Carbonate Aquifer Void Evolution (CAVE) code previously published by Teutsch (1993), Sauter (1993), Clemens et al. (1996), Bauer et al. (2000, 2003), Liedl et al. (2003), and Birk (2002). Thus, this article is focused on the new algorithm developed for CFPM2 (Shoemaker et al. 2008). Although Jeannin (2001) modified the Darcy-Weisbach equation and described turbulent flow to a spring with the one-dimensional Equation 2, sffiffiffiffiffiffiffiffiffiffiffiffi 
H  Qturbulent ¼ k#s  ð2Þ
l  where k#s is the conductance [units of L3T21] of the spring and
H/
l is the absolute value of the gradient [dimensionless], the use of Equation 2 would require the complete reformulation of a solution of a different partial differential equation for flow, and MODFLOW could not be used readily. Halford (2000) showed that turbulent and laminar flow in well bores can be solved in terms of the Darcy equation if hydraulic conductivity is made a nonlinear function of the laminar and turbulent friction factors for pipes allowing only minor modifications to MODFLOW. Our approach for simulation of laminar and turbulent flow in porous media for CFPM2 is similar to the approach taken in Halford (2000).

Specific velocity under turbulent conditions diverges from Darcy’s law proportionately to the square root of specific velocity. The modified form of Darcy’s law that is incorporated in CFPM2 is as follows: rffiffiffiffiffi vC
h ; v > vC ð3Þ v ¼ KLAM v
L where KLAM is the laminar hydraulic conductivity, vC is the specific velocity at the onset of turbulent flow, and
h/
L is the hydraulic gradient. Specific velocity is solved iteratively because specific velocity also is on the right-hand side of Equation 3. Turbulent flow is simulated in CFPM2 by computing a turbulent hydraulic conductivity (KTURB) as a nonlinear function of the Re, when Re is greater than NRe. Turbulent loss is defined with Reynolds number because differences in mean pore diameter also affect specific or mean velocity. The equation relating KTURB and Re is as follows: rffiffiffiffiffiffiffiffi NRe KTURB ¼ KLAM ð4Þ ; Re > NRe Re The KTURB equals KLAM when Re is less than or equal to NRe (Shoemaker et al. 2008). A power function approach similar to Equation 4 has been applied previously to simulate turbulent flow in well completions and boreholes (Halford 2000). Hydraulic conductivity is reduced by the square root of the laminar friction function divided by the turbulent friction function in the turbulent flow module of MODFLOW-96 (Halford 2000). Laminar and turbulent pipe-flow friction factors are defined with the Darcy-Weisbach and Sacham formulas, respectively (Geankoplis 1993). Using the simpler approach in CFPM2 (Equation 4) appeared adequate for simulating turbulent head losses because pipe roughness minimally affects friction factors (Halford 2000). Additionally, empirical equations for pipe friction factors are functions of the Re. The power function in CFPM2 (Equation 4) is solved iteratively because Re is a function of specific velocity. Spreadsheet solutions of Equation 3 require solving for KTURB and specific velocity in multiple cells or iterating in a macro. Specific velocity, heads, KTURB, and Re will converge in CFPM2 if an iterative solver such as the strongly implicit procedure or preconditioned conjugate gradient is used. Turbulent head loss computations are triggered in CFPM2 when head differences between cells exceed critical thresholds rather than explicitly calculating Re and comparing to NRe. Critical head difference is calculated in CFPM2 as:
hCRIT ¼

NRe l
L qKLAM d

ð5Þ

where
L is the row or column spacing. Critical head difference is tracked because all other variables in the Re calculation (Equation 1) and cell-to-cell conductances are invariant while flow is laminar. Critical head differences are calculated between rows and columns of cells in conduit flow layers (Shoemaker et al. 2008). E.L. Kuniansky et al. GROUND WATER 46, no. 5: 768–771

769

40

The RMS error for all 24 observations of specific velocity was 0.6 m/d. Flow in cube 6 ceased being laminar when the specific velocity exceeded 11 m/d. The maximum difference between KTURB and KLAM occurred along horizontal axis 1 because KLAM was greatest along this axis, which created higher specific velocities (Figure 1b). Critical hydraulic head differences differed for each direction because KLAM differed between axes, as predicted with Equation 3. The verification of the simple empirical power function in Equation 4 applied in CFPM2 is with laboratory data from porous media consisting of carbonate rock samples with touching vug porosity and pore diameters of 0.01 m and at fairly low Reynolds numbers.

(a) Cube 6 Vertical Axis

30

CRITICAL VELOCITY MEASURED SIMULATED

20

10

Klam = 48 m/d NRe = 1.44 RMS = 0.70 m/d

SPECIFIC VELOCITY, IN METERS PER DAY

0 40

(b) Cube 6 Horizontal Axis 1

Guidance on Assignment of Parameters for Mode 2

30

20

10

Klam = 61 m/d NRe = 1.44 RMS = 0.35 m/d

0 40

(c) Cube 6 Horizontal Axis 2

30

20

10

Klam = 36 m/d NRe = 1.44 RMS = 0.57 m/d

0 0.00

0.25

0.50

0.75

1.00

∆H/∆L, DIMENSIONLESS

Figure 1 a-c. Measured and simulated (CFPM2) relations between gradient DH/DL where DL is 0.2 m, and specific velocity for directions of cube 6. (Measured data from DiFrenna et al. 2007, figure 6).

Verification of Power Function in CFPM2 The measured experimental data from a 0.2-m cube (DiFrenna et al. 2007, figure 6) and results simulated using the power function applied in the CFPM2 agree closely (Figure 1). A laminar hydraulic conductivity (KLAM) was estimated for each axis, and a single NRe was estimated for cube 6 first using a spreadsheet (four parameters estimated with 24 observations). The root mean square (RMS) error between simulated and measured specific velocities was minimized, using least squares regression, to estimate the three KLAM and NRe. The KLAM estimates were 48, 61, and 36 m/d for the vertical, horizontal 1, and horizontal 2 axes, respectively. The NRe of 1.4 was estimated for all directions of cube 6. 770

E.L. Kuniansky et al. GROUND WATER 46, no. 5: 768–771

Simulating turbulent flow with the CFP requires more data than traditional MODFLOW simulations. CFPM1 is an update of subroutines of the CAVE code previously published by Teutsch (1993), Sauter (1993), Birk (2002), Clemens et al. (1996), Liedl et al. (2003), and Bauer et al. (2000, 2003) to MODFLOW-2005. Guidance on the use of CFPM1 can be found in the previous references and in Shoemaker et al. (2008). For CFPM2, additional data requirements include the temperature of ground water in the preferential flow layers, mean void diameters, and lower and upper critical Reynolds numbers (NRe) for all the preferential flow layers and the laminar horizontal hydraulic conductivities, KLAM. Assigning values for these data requirements can be challenging. Additionally, one cannot use a quasithree-dimensional modeling approach for CFPM2 preferential flow layers because the top and bottom of the preferential flow layer must be assigned as part of the MODFLOW-2005 discretization file input. Top and bottom elevations for preferential flow layers can be determined using geologic mapping, well lithology data, and borehole and/or surface geophysics. Again, the temperature of ground water should be relatively easy to quantify using water temperature probes. Mean void diameters could be assigned based on visual inspection of geological cores or images of the borehole. The upper and lower critical Reynolds numbers and laminar horizontal hydraulic conductivity, KLAM, are the most difficult properties to assign for CFPM2 simulations. Values for upper and lower critical Reynolds numbers and laminar horizontal hydraulic conductivities depend on the geology of each preferential flow zone. When borehole flow meter logging is used in conjunction with traditional aquifer tests, an estimate of laminar hydraulic conductivity can be obtained for the preferential flow layers (Paillet 1998; Reese and Cunningham 2000; Cunningham et al. 2006). Based on the experimental data published by DiFrenna et al. (2007), if the interconnected vugs are small (0.01 m or less), the upper and lower critical Reynolds numbers are close to that of porous media (1 to 50). For larger interconnected vugs, the critical Reynolds numbers are probably much less than 2000, the lower critical Reynolds number for pipes (Vennard and Street 1975).

Conclusions Turbulent flow in aquifers can be simulated with a simple power function (Equation 4). The transition from laminar to turbulent (non-Darcian) flow is simulated well with this simple approach. Experimental verification was limited to Reynolds numbers less than three times the critical Reynolds number. The approach should extend to greater Reynolds numbers because specific velocity is proportional to the square root of hydraulic gradient where flow is fully turbulent. The simple power function in the CFPM2 (Equation 4) can be applied more easily than complex pipe-flow equations. Pipe-flow equations require an additional roughness term (Halford 2000). Additionally, these equations can be unstable when NRe approaches 1. The relation between gradient and specific velocity is solved equally well using either the CFPM2 approach or pipeflow equations when NRe greater than 10. For simulation of larger submerged conduit systems, CFPM1 should be applied. However, the use of CFPM1 requires that the user has maps or the network of caverns, ground water temperature, wall roughness, tortuosity, and critical Reynolds numbers for each simulated pipe segment. MODFLOW-2005 with the CFP can be downloaded from a link on the MODFLOW-2005 Web site page, http://water.usgs.gov/software/ground_water.html.

Acknowledgments The authors express their appreciation to Andrew Long and Randall Laczniak of the USGS for their internal review prior to submission to the journal and Steffen Birk, Gareth Davies, and an anonymous reviewer of the journal who helped improve the clarity and content of the paper.

References Bauer, S., R. Liedl, and M. Sauter. 2003. Modeling of karst aquifer genesis: Influence of exchange flow. Water Resources Research 39, no. 10: 1285. Bauer, S., R. Liedl, and M. Sauter. 2000. Modelling of karst development considering conduit-matrix exchange flow. In Calibration and Reliability in Groundwater Modelling: Coping with Uncertainty—Proceedings of the ModelCARE’99 Conference, ed. F. Stauffer, W. Kinzelbach, K. Kovar, and E. Hoehn, 10–15. International Association of Hydrological Sciences Publication 265. Wallingford, Oxfordshire, UK: International Association of Hydrological Sciences. Bear, J. 1979. Hydraulics of Groundwater. New York: McGrawHill Inc. Birk, S. 2002. Characterization of karst systems by simulating aquifer genesis and spring responses: Model development and application to gypsum karst. Vol. 60 of Tu¨binger Geowissenschaftliche Arbeiten: Tu¨bingen, Germany, Reihe C. Institut und Museum fu¨r Geologie und Pala¨ontologie der Universita¨t Tu¨bingen. Clemens, T., D. Hckinghaus, M. Sauter, R. Liedl, and G. Teutsch. 1996. A combined continuum and discrete network reactive transport model for the simulation of karst development. In Calibration and Reliability in Groundwater

Modeling, ed. K. Kovar and P. van der Heijde, 309–318. International Association of Hydrological Sciences Publication 237. Wallingford, Oxfordshire, UK: International Association of Hydrological Sciences. Cunningham, K.J., M.A. Wacker, E. Robinson, J.F. Dixon, and G.L. Wingard. 2006. A cyclostratigraphic and borehole geophysical approach to development of the three-dimensional conceptual hydrogeologic model of the karstic Biscayne aquifer, southeastern Florida. USGS Scientific Investigations Report 2005-5235. Reston, Virginia: USGS. Darcy, H. 1856. Les Fontances publiques de la ville de Dijon. Paris: Victor Dalmont. [English translation by Patricia Bobeck. 2004, Kendall/Hunt Publishing.] DiFrenna, V.J., R.M. Price, and M.R. Savabi. 2007. Identification of a hydrodynamic threshold in karst rocks from the Biscayne Aquifer, south Florida, USA. Hydrogeology Journal, doi 10.1007/s10040-007-0219-4. Freeze, R.A., and Cherry, J.A. 1979. Groundwater. Englewood Cliffs, New Jersey: Prentice-Hall. Geankoplis, C.J. 1993. Transport Processes and Unit Operations, 3rd ed. Englewood Cliffs, New Jersey: Prentice-Hall. Halford, K.J. 2000. Simulation and interpretation of borehole flowmeter results under laminar and turbulent flow conditions. In Proceedings of the Seventh International Symposium on Logging for Minerals and Geotechnical Applications, Golden, Colorado, October 24–26, 2000, 157–168. Houston, Texas: The Minerals and Geotechnical Logging Society. Harbaugh, A.W. 2005. MODFLOW-2005, the U.S. Geological Survey modular ground-water model—The ground-water flow process. USGS Techniques and Methods 6-A16. Reston, Virginia: USGS. Jeannin, P.-Y. 2001. Modeling flow in the phreatic and epiphreatic karst conduits in the Holloch cave (Muotatal, Switzerland). Water Resources Research 37, no. 2: 191–200. Liedl, R., M. Sauter, D. Hu¨ckinghaus, T. Clemens, and G. Teutsch. 2003. Simulation of the development of karst aquifers using a coupled continuum pipe flow model. Water Resources Research 39, no. 1: 1057, doi:10.1029/ 2001WR001206. Paillet, F.L. 1998. Flow modeling and permeability estimation using borehole flow logs in heterogeneous fractured formations. Water Resources Research 34, no. 5: 997–1010. Reese, R.S., and K.J. Cunningham. 2000. Hydrogeology of the Gray Limestone Aquifer in southern Florida. USGS WaterResources Investigations Report 99-4213. Restn, Virginia: USGS. Sauter, M. 1993. Double porosity models in karstified limestone aquifers: Field validation and data provision. In Hydrogeological Process in Karst Terranes, 261–279. Wallingford, Oxfordshire, UK: International Association of Hydrological Sciences 207. Schneebeli, G. 1955. Experiences sur la limite de validite de la loi de Darcy et l’apparition de la turbulence dans un ecoulement de filtration. La Houille Blanche 10, no. 2: 141–149. Shoemaker, W.B., E.L. Kuniansky, S. Birk, S. Bauer, and E.D. Swain. 2008. Documentation of a conduit flow process (CFP) for MODFLOW-2005. USGS Techniques and Methods 6-A24. Reston, Virginia: USGS. Teutsch, G. 1993. An extended double-porosity concept as a practical modeling approach for a karstified terraine. In Hydrogeological Processes in Karst Terraines, Proceedings of the Antalya Symposium and Field Seminar, October 1990, ed. J. Gultekin and Back, 281–292. International Association of Hydrological Sciences Publication 207. Wallingford, Oxfordshire, UK: International Association of Hydrological Sciences. Vennard, J.K., and R.L. Street. 1975. Elementary Fluid Mechanics, 5th ed. New York: John Wiley.

E.L. Kuniansky et al. GROUND WATER 46, no. 5: 768–771

771