Bench scale laboratory tests to analyze non-linear flow in fractured media

. The knowledge of ﬂow phenomena in fractured rocks is very important for groundwater resources management in hydrogeological engineering. A critical emerging issue for fractured aquifers is the validity of the Darcian-type “local cubic law”, which assumes a linear relationship between ﬂow rate and pressure gradient to accurately describe ﬂow patterns. Experimental data obtained under controlled conditions such as in a laboratory increase our understanding of the fundamental physics of fracture ﬂow and allow us to investigate the presence of non-linear ﬂow inside fractures that generates a substantial deviation from Darcy’s law. In this study the presence of non-linear ﬂow in a fractured rock formation has been analyzed at bench scale in laboratory tests. The effects of non-linearity in ﬂow have been investigated by analyzing hydraulic tests on an artiﬁcially created fractured rock sample of parallelepiped (0 . 60 × 0 . 40 × 0 . 8 m) shape. The volumes of water passing through different paths across the fractured sample for various hydraulic head differences have been measured, and the results of the experiments have been reported as speciﬁc ﬂow rate vs. head gradient. The experimental results closely match the Forchheimer equation and describe a strong inertial regime. The results of the test have been interpreted by means of numerical simulations. For each pair of ports, several steady-state simulations have been carried out varying the hydraulic head difference between the inlet and outlet ports. The estimated linear and non-linear Forchheimer coefﬁcients have been correlated to each other and respectively to the tortuosity of the ﬂow paths. A correlation among the linear and non-linear Forchheimer coefﬁcients is evident. Moreover, a tortuosity factor that in-ﬂuences ﬂow dynamics has been determined.


Introduction
The purpose of this paper is to experimentally investigate the behavior of high velocity flow regimes in a fractured network at bench scale.
High velocity flow dynamics can have significant impact in diverse fields such as radioactive waste disposal, geothermal engineering, environmental risk assessment and remediation, reservoir engineering, and groundwater hydrology (Cherubini et al., 2010).
In most studies examining hydrodynamic processes in saturated porous and fractured media, it is assumed that flow is described by Darcy's law, which expresses a linear relationship between pressure gradient and flow rate.Darcy's law has been demonstrated to be valid at low flow regimes (Re ≪ 1).For Re > 1 a non-linear flow behavior is likely to occur.Ing and Xiaoyan (2002) have shown how non-Darcian flow has a significant impact on consolidation rate in geotechnics.Basak and Rajagopalan (1982) have demonstrated that seawater intrusion length increases when flow deviates from Darcian linearity.
Exposure and risk assessment of chemical pollution based on the applicability of Darcy's law are sometimes inappropriate when applied to fractured aquifers because of the Published by Copernicus Publications on behalf of the European Geosciences Union.
possible occurrence of non-laminar flow regimes.In fact, considerable evidence exists to refute their reliability in assessing solute-migration rates and for determining downgradient concentrations (Field, 1997).In solute transport, non-Darcian flow might give rise to tailing in contaminant breakthrough curves that show a non-Fickian behavior (Boutt et al., 2006;Cardenas et al., 2007).
A non-Fickian behavior is also valid in case of highconcentration brine transport, where the assumption is that the linear Darcy's law holds as shown in Watson et al. (2002).On the other hand, Tenchine and Gouze (2005) carried out density driven flow simulations in a rough-walled natural fracture extracted from a limestone quarry and observed measurable non-linearity between the velocity growth rate and the velocity, indicating that the Navier-Stokes formulation must be used.
The mathematical representation of fluid dynamics in fractured rock aquifers is of a great concern for environmental and petroleum engineering and in geological sciences.
The local cubic law adapts Darcy's law to flow through fractures under the assumption of ideal fractures with flat, smooth and parallel walls and infinite lengths, together with laminar flow, incompressible fluid and confined aquifer configuration.
Real rock fractures instead are characterized by rough walls and variable surfaces, geometry and apertures.The presence of asperities and obstructions or sharp changes in fracture profile is the reason for microscopic inertial phenomena that cause an extra macroscopic pressure loss which deviates flow from linearity.
Tortuosity, which usually characterizes the ratio of the effective path length connecting two locations in porous media to the geometric distance (Tenchine and Gouze, 2005), has been found to significantly affect fluid flow in fractured media under certain conditions (Tsang, 1984;Wang and Narasimhan, 1985).Yeo and Ge (2005) identified a criterion parameter function of the roughness and tortuosity for the applicability of the Reynolds equation to fluid flow in rock fractures.Louis (1968) empirically defined five steady-state flow regimes, from smooth laminar to fully developed turbulent, depending on various degrees of relative roughness, fracture aperture, hydraulic head gradient in the fracture plane and kinematic viscosity.
Fully developed turbulent flow conditions are more likely to occur when velocities are high, i.e. in karstic conduits, whereas, in presence of lower velocities, non-linear laminar flow happens when inertial effects become important, e.g.due to roughness, restrictions, presence of contact points between the fracture walls and non-rectilinear profile.
Different from conduits, the flow regime prevailing between rough-walled surfaces does not immediately switch from Darcy flow to turbulent flow as velocity is increased or vice versa (Kohl et al., 1997), but transitional regimes can be distinguished.
On the basis of a recent classification by Skjetne (1995), three different principal flow regimes can be distinguished: (1) linear laminar flow, (2) non-linear rough laminar flow (weak or strong inertia) and (3) turbulent flow.
A fracture can be envisioned as two rough surfaces in contact.Cross sectional solid areas representing asperities in contact are similar to the grains of porous media.It is therefore possible to apply the general equations describing flow and transport in porous as well as in fractured media.The void space of the fracture is interconnected in one surface, so the fracture is best described as a two-dimensional porous medium (National Research Council, 1996).
In the literature different laws are reported that account for the non-linear relation between velocity and pressure gradient.
The weak inertia equation is a cubic extension of Darcy's law and describes pressure loss versus flow rate for low flow rates: where ) is the velocity, and γ (L) is called the weak inertia factor.This law was first shown numerically (Barrére, 1990) and then derived theoretically for homogeneous isotropic media (Mei and Auriault, 1991).
In case of higher Reynolds numbers (Re ≫ 1), the pressure losses pass from a weak inertial to a strong inertial regime, described by the Forchheimer equation (Forchheimer, 1901), given by: where β is called the inertial resistance (coefficient), or non-Darcy coefficient.
Forchheimer proposes a further expression in which the pressure gradient can be described by a third order polynomial, whereas Hassanizadeh and Gray (1987) propose an extension that includes the effects of a transient flow regime.
Another commonly used non-Darcian flow equation is the Izbash or power-law equation (Izbash, 1931): where λ and n are two constant coefficients and 1 ≤ n ≤ 2. The Forchheimer equation can be considered as an expansion of the linear Darcy's law because it includes both the viscous and inertial effects.
On the other hand, power law functions such as the Izbash equations are more appropriate to model post-linear non-Darcian flows that might be caused predominately by turbulent effects (Wen et al., 2008).
Experimental evidence has shown that both the Forchheimer and the Izbash equations are equally possible, depending on field-specific conditions (Wen et al., 2006).
A general Darcian-like relationship can be used (Chin et al., 2009) to describe all the mentioned flow regimes and to account for non-linearities in the relationship between hydraulic head gradient and flow velocity: where K eff (L T −1 ) is the effective hydraulic conductivity, and h = p/ρg (L) is the total hydraulic head.For instance, according to Forchheimer's law, the effective hydraulic conductivity can be written as (Cakmak, 2009): where a (T L −1 ) and b (T 2 L −2 ) are the linear and inertial coefficients, respectively, in terms of hydraulic head and can be expressed according to the previous expression as: In the same way the effective fracture transmissivity for a discontinuity can be defined, and the relationship between specific flow rate q (L 2 T −1 ) and hydraulic head gradient can be written as: where a f and b f are related to a and b: where w represents fracture aperture.Flow regimes and non-linear behavior of fluid flow through fractures have been investigated empirically (Lomize, 1951;Louis, 1969), experimentally (Witherspoon et al., 1980;Qian et al., 2005;Chin et al., 2009), andnumerically (Zimmerman et al., 2004;Kolditz, 2001;Brush and Thomson, 2003).
Zimmerman et al. ( 2004) studied non-linear flow regimes both with experimental and numerical simulations of Navier-Stokes equations and suggested the critical Reynolds number of 10 for practical purposes.2.44 Hydraulic conductivity (m s −1 ) 1.63e-8 Witherspoon et al. (1980) conducted experiments in rock fractures to study the hydraulic behavior.They demonstrated that if fractures are rough or flow velocities are high, Darcy's law is not applicable.
Chin et al. ( 2009) investigated the variation of effective hydraulic conductivity as a function of specific discharge in several 0.2 m and 0.3 m cubes of Key Largo Limestone.The experimental results closely match the Forchheimer equation.Qian et al. (2005) carried out laboratory experiments to study groundwater flow in a single fracture under different conditions of fracture aperture, surface roughness, and water pressure.The experimental results show that the average flow velocity (V ) could be approximated by an empirical exponential function of the hydraulic gradient (I ), which varies in the range of 0.003-0.02.The power index of the exponential function is close to 0.5.Such a V -I relationship indicates a non-Darcian turbulent flow, even though the Reynolds number is relatively low.Javadi et al. (2010) performed both laminar and turbulent flow simulations for a wide range of flow rates in an artificial three-dimensional fracture.They developed a new geometrical model for non-linear fluid flow through rough fractures, which suggests a polynomial expression, like the Forchheimer law, to describe the dependence of pressure drop on flow rate.Finally, this model has been evaluated with experimental results of a fracture with different geometries.A good accuracy was found between the proposed model and turbulent flow simulation results.Brush and Thomson (2003) carried out simulations of fluid flow through single rough-walled fractures with various apertures using Navier-Stokes, Stokes, and local cubic law simulations.They observed that several Navier-Stokes velocity profiles have flatter peaks or noses that indicate the formation of an inertial core between the walls.They demonstrated that inertial forces can significantly influence the internal flow field within a fracture, so that the forces driving the flow field are reduced, and the overall flow rate is decreased.As the mean aperture increases, the effect of surface roughness diminishes (Boutt et al., 2006).
In the present paper, several hydraulic tests have been carried out in a fractured rock sample, and a non-linear behavior of the flow regime has been ascertained.The study is aimed at determining the relationship between the average velocity and the hydraulic gradient for each investigated fracture-controlled flow pathway, and will serve as the first  2 Experimental setup

Characterization and preparation of block sample
The experiments have been performed on a limestone block with parallelepiped shape (0.6 × 0.4 × 0.08 m 3 ) recovered from the "Calcare di Altamura" formation, which is located in the Apulia region of southeastern Italy.
In Table 1 are shown the bulk hydraulic parameters of the limestone block.The fracture network has been made artificially through 5 kg mallet blows.The fissured system and fracture apertures on the block surfaces have been recorded with a high resolution digital camera.Subsequently, the images have been scaled and rectified using Perspective Rectifier (www.rectifiersoft.com)(Fig. 1).Profiles of discontinuities and aperture measurements have been extracted from the recorded images using edge function with the "canny" filter from the built-in Scilab Image Processing Toolbox (www.scilab.org).For each discontinuity the median profile and the aperture distribution have been determined (Table 2).
The surface of the block sample has been sealed with transparent epoxy resin (Leven et al., 2004).A hole of 1 cm diameter has been opened for each discontinuity in correspondence to the boundary of the block.
The spatial position of opened ports (Fig. 2) and the illustration of construction details (Fig. 3) have been reported.

Materials and methods
In Fig. 4 the diagram of the experimental set up is shown.The sealed block sample is connected with a hydraulic circuit.Water moves from the upstream to the downstream tank and returns to the upstream tank by means of a transfer pump.A flow cell is connected to the outlet port.The sealed block and the tubes of hydraulic circuits are completely saturated.Initially, the valves "a" and "b" are closed, and the hydrostatic head in the flow cell is equal to h 0 .The ultrasonic flow velocimeter measures the snapshot flow rates that enter the sealed block.The experiment begins with the opening of the valve "a", which is reclosed when the hydraulic head in the flow cell is equal to h 1 .Finally, the hydraulic head in the flow cell is reported to h 0 through the opening of the valve "b".
The average flow rate through the sealed block can be estimated by means of the volumetric method: where S 1 (L 2 ) represents the cross-sectional area of the flow cell.During the experiments, these values are compared with the snapshot flow rates measured by the ultrasonic velocimeter in order to check the absence of hydraulic losses due to obstructions and leaks in the hydraulic system.The experiment is repeated changing the hydraulic head h c (control head) of the upstream tank, and for each configuration of inlet-outlet ports.For different values of h c , the time, t = (t 1 − t 0 ), required to fill the flow cell from h 0 to h 1 has been registered.
The storage property of the upstream tank (S 2 ) is much higher than that of the downstream flow cell (S 2 ≫ S 1 ), therefore during the experiment the upstream hydraulic head can be considered constant.Given that the sealed block sample and the hydraulic circuit are very rigid, their compressibility can be neglected.
On the basis of these assumptions the drainage process is governed by the following equation:  where h (L) is the hydraulic head of the downstream flow cell, h c (L) is the hydraulic head of the upstream tank, and Ŵ ( h) (L 2 T −1 ) is the hydraulic conductance term (Harbough, 2005) representative of both hydraulic circuit and the active fracture network configuration.Hydraulic losses within the single hydraulic circuit can be expressed according to Chezy law as: where C is a characteristic coefficient related to the roughness, section and length of the tubes of the hydraulic circuit.Whereas, only for the sealed block, the h-Q relationship can be expressed through the following polynomial expression: where A and B are the linear and non-linear hydraulic loss coefficients, respectively, and are related to the roughness, aperture, length and shape of fractures.Combining Eqs. ( 11) and ( 12), the conductance term representative of the whole hydraulic system assumes the following expression: Substituting Eq. ( 13) into Eq.( 10) and integrating the latter from t = t 0 to t = t 1 with the initial condition h = h 0 , the following equation is obtained: Then, fitting the experimental relationship between the time, t = (t 1 − t 0 ), and the hydraulic head of the downstream

Experimental results
Several experiments have been carried out for each in-out port configuration.The control head h c varies in the range of 0.17-1.37m, and the average flow rates observed are in the range of 3.08 × 10 −7 -2.99 × 10 −5 m 3 s −1 .All the executed experiments show a non-linear Q-h relationship that can be well described by Eq. ( 12). Figure 5 shows the fitting method described in the previous section to estimate the linear (A) and non-linear (B) terms.The double entry Table 3 shows the estimates of A and B for each pair of ports.Figure 6 shows flow rate Q versus hydraulic head difference dh for all the experiments.If the experiments were carried out on a single fracture, a f and b f coefficients of Forchheimer Eq. ( 5) could be derived in an analytical way from Eq. ( 12).In the present case, in order to obtain Forchheimer terms for each path a numerical model has been implemented.
Starting from fracture profiles (Fig. 1), the threedimensional geometry of the fissured system has been reconstructed using the GMSH pre-processing tool (Geuzaine and Remacle, 2009).The geometry has been imported in Comsol Multiphysics ® v 4.0a (Comsol AB, 2010) using STL exchange file format (Vinciguerra and Bernabè, 2009).Furthermore, the geometry of the port holes has been modeled.
COMSOL uses the finite element scheme to solve generic partial differential equations.In particular, "Weak Form Boundary PDE Interfaces" included in "Mathematics Module" has been used to solve the following continuity equation: Figure 7a shows the mesh of the finite element model used for numerical simulations.a f and b f are constant for the whole domain and represent the equivalent parameters for linear and non-linear hydraulic losses.The aim of the numerical model is to estimate a f and b f for each port configuration.By means of Eq. ( 9), an estimate of flow rate Q obs for steady-state conditions can be obtained in correspondence to a hydraulic head difference  imposed between the inlet and outlet ports.In a similar way by means of numerical simulations, a relationship between the hydraulic head difference imposed between inlet and outlet ports and the simulated flow rate Q sim in correspondence to outlet ports can be obtained by varying a f and b f .The idea is to find the equivalent parameters a f and b f for each path that minimizes the difference between Q obs , evaluated through Eq. ( 9), and Q sim , evaluated by means of numerical model for each imposed hydraulic difference.The double entry Table 4 shows a f and b f estimated for each pair of ports. 1 -2.96e3 6.34e3 2.64e3 5.29e3 3.11e3 3.02e4 2 1.94e7 -8.04e3 3.83e3 6.30e3 3.88e3 2.99e4 3 5.77e7 1.02e8 -5.86e3 8.43e3 6.97e3 3.00e4 4 9.32e6 2.91e7 6.33e7 -6.10e3 2.97e3 2.06e4 5 5.69e7 8.34e7 1.42e8 7.15e7 -1.33e4 2.62e4 6 1.66e7 3.39e7 8.92e7 1.51e7 3.80e8 -3.91e4

Discussion
In order to analyze the experimental results, two dimensionless numbers have been evaluated: the Reynolds number and the Forchheimer number.
Reynolds number is defined as the ratio of inertial forces to viscous forces: where v represents the average velocity evaluated on the active path, and D h represents the characteristic dimension.For a fracture having small aperture with respect to its height, the characteristic dimension radius is equal to D h = w/2.Assuming that the average specific discharge is equal to q = v • w, the Reynolds number becomes: The Forchheimer number is the ratio of non-linear to linear pressure losses: Reynolds number versus Forchheimer number for all paths.According to Eq. ( 8) the Forchheimer number can be reformulated as: According to Zeng and Grigg (2006), the Forchheimer number is recommended as a criterion for identifying non-Darcy flow because it has the advantage of clear meaning.It equals the ratio of pressure drop caused by liquid-solid interactions to that by viscous resistance, and it is directly related to the non-Darcy effect.Inertial effects dominate over viscous effects at the critical Forchheimer number (Fo > 1) (Ruth and Ma, 1992).The Reynolds number is instead a dimensionless number that indicates when microscopic inertial effects become important.It is inappropriate on the macroscopic level, because microscopic inertial effects do not directly lead to macroscopic inertial effects.In fact, a high microscopic Reynolds number does not necessarily imply non-Darcian flow.Instead, Fo indicates precisely the onset of non-Darcian flow (Ruth and Ma, 1992); it accounts for both velocity (v) and structure of the medium because β is structure dependent.The term β inherently contains information on the tortuosity of the flow paths that leads to changes in the microscopic inertial terms.In fact, if the structure of the medium is such that microscopic inertial effects are rare, then β will be small, and Fo will remain small until v (i.e.Re) is large.Instead, both β and Fo will be large if the structure of the medium is such that microscopic inertial effects can be expected.
In even if the paths (1-4 and 4-6) reach relatively high values of Re, they present values of Fo lower than unity.
The slope of the straight lines represents the ratio between Re and Fo that is equal to the ratio of their respective characteristic dimensions: This dimensionless group ζ is characteristic of the flow path.
Relatively high values of this parameter correspond to a more linear flow behavior, because inertial effects dominate viscous ones at higher Reynolds numbers.Therefore it permits us to distinguish a different behavior of the experiments carried out with varying configurations of ports.
In particular, in Fig. 8 it is possible to distinguish a set with shallower slopes.This set has in common the outlet port 7.In Fig. 9 the shape of the fracture in correspondence to port 7 is shown.The particular shape of this fracture gives rise to a greater contact between fracture surfaces compared to the others.In fact, the path that contains this fracture presents a very high hydraulic loss.
For each path the equivalent aperture w eq has been estimated from the linear term assuming the cubic law is valid: This term has been compared with the average measured aperture of each path w (Fig. 9).Though w eq underestimates w, they are of the same order of magnitude.
A power law has been observed between the terms a and b (Fig. 11): This correlation between inertial and viscous coefficients is customarily used in petroleum production engineering in order to predict high velocity well performance.Geertsma (1974) and Skjetne et al. (1999) found similar relationships for high velocity flow in porous media and fractures.
Zimmerman and Bodvarsson (1996) analyzed flow in twodimensional rough-walled rock fractures and found a factor equal to the ratio of the cubes of w eq and w that reflects the tortuosity induced into the streamlines by the obstacles.
In the case of a single rock fracture, this factor depends only on the area of the asperity region.However, in the case study it depends on several parameters such as roughness, the areas where the rock faces are in contact with each other, fracture intersection, position and shape of the inlet and outlet ports: This parameter measures how much each path deviates from the parallel plate model.The pressure drops depend by the morphology of the fracture wall surfaces and on the tortuosity of the flow paths.Significant head losses may be envisaged to occur adjacent to sharp corners of a fracture where sudden change of aperture occurs (Javadi et al., 2010).
Figures 12 and 13 show τ -b and τ -a relationships, respectively.τ results correlated with b and a by means of a power function.

Conclusions
In this paper non-linear fluid flow through rock fractures was studied by means of laboratory tests and numerical modeling.
The Forchheimer equation has proven highly effective in explaining the relationship between flow rate and pressure drop, which depicts a strong inertial regime where the viscous and inertial pressure drops are controlled by v and v 2 term, respectively.
The equivalent linear and non-linear terms of Forchheimer's law have been estimated by numerical modeling.
The equivalent aperture of each flow path is determined assuming the cubic law is valid.Though it underestimates the mean measured aperture, it keeps its same order of magnitude.
A tortuosity factor τ has been determined as the ratio of the cube of the equivalent aperture and the cube of the mean measured aperture.This factor is a measure of the deviation of each flow path from the parallel plate model.In other words it measures the effects of different factors such as roughness, the contact area between fracture surfaces, fracture intersections, and the position and the form of the inlet and outlet ports.
A power law between the Forchheimer terms and τ has been detected.In complex fracture networks the tortuosity factor plays an important role in fluid flow dynamics.
The experimental results showed that the dependence of hydraulic conductivity on specific discharge cannot be neglected in fractured media.For instance, during pumping tests a linear flow model can cause errors in the determination of transmissivity in fractured rock aquifers, because much of the data collected can be non-linear due to flow occurring in transition between linear and fully turbulent flow.
On the other hand, potential errors induced by a non-linear flow model in constant pressure tests have also been recognized in the engineering literature (Louis and Maini, 1970).Elsworth and Doe (1986) used mathematical modeling of packer tests in fractured rock to show that calculation of transmissivity using non-Darcian constant head data can lead to underestimation errors as much as an order of magnitude.
In our study the effective hydraulic transmissivity proves to be less than 46.59 % (average value) of the Darcian (linear) flow hydraulic transmissivity.In particular, in correspondence to path 3-4 the variation reaches a maximum value equal to 59.38 %.
In pumping tests multiple pressure steps (i.e. higher flow rates resulting in a much greater differential pressure) should be used for a more accurate identification of the Darcian range and the quantification of the linear to non-linear flow relations, resulting in better transmissivity estimates as a non-linear function of the gradient.
This concept has to be taken into account in cases of anthropogenic stresses in the aquifer that might give rise to high hydraulic gradients.

kFig. 2 .
Fig. 2. Localization of the ports on the top surface of the block; (a) top fracture traces, (b) bottom fracture traces, and (c) ports.

Figure 4 .
Figure 4. Schematic diagram of the experimental setup (picture is not to scale).516 Fig. 4. Schematic diagram of the experimental setup (picture is not to scale).

Fig. 5 .
Fig. 5. Experimental results obtained for 1-2 configuration ports.(a) Control head H c versus time.(b) Average flow rate Q versus difference head evaluated as dh = H c + (h 0 + h 1 )/2.(c) Difference head versus conductance term evaluated as Eq.(13).(d) Average flow rate versus resistance term evaluated as the inverse of conductance.Dots represent the experimental values, dashed line represent the fitting of experimental values, and marked line represent the functions without the effect of circuit (C = 0).

523Figure 6 Fig. 6 .
Figure 6 Flow rate Q versus hydraulic head difference dh for all the experiments.Dots represents the experimental values, Fig. 7. (a) Finite element mesh of the numerical model.(b) A result of a simulation, color scale indicates hydraulic head, black arrows represent specific discharge.

Fig. 9 .
Fig. 9. Fractured block.The fracture in correspondence to the hole 7 is shown in the lower left side.

Table 1 .
Properties of the limestone block.

Table 2 .
Characteristic parameters of the discontinuities.

Table 3 .
Double entry table.The upper and lower triangular matrix represents the linear coefficient (A) and the non-linear coefficient (B), respectively, obtained for each pair of the ports.

Table 4 .
Double entry table.The upper and lower triangular matrix represents the linear coefficient (a f ) and the non-linear coefficient (b f ), respectively, obtained for each pair of the ports.