Evidence of non-Darcy flow and non-Fickian transport in fractured media at laboratory scale

During a risk assessment procedure as well as when dealing with cleanup and monitoring strategies, accurate predictions of solute propagation in fractured rocks are of particular importance when assessing exposure pathways through which contaminants reach receptors. Experimental data obtained under controlled conditions such as in a laboratory allow to increase the understanding of the fundamental physics of fluid flow and solute transport in fractures. In this study, laboratory hydraulic and tracer tests have been carried out on an artificially created fractured rock sample. The tests regard the analysis of the hydraulic loss and the measurement of breakthrough curves for saline tracer pulse inside a rock sample of parallelepiped shape (0.60× 0.40× 0.08 m). The convolution theory has been applied in order to remove the effect of the acquisition apparatus on tracer experiments. The experimental results have shown evidence of a nonDarcy relationship between flow rate and hydraulic loss that is best described by Forchheimer’s law. Furthermore, in the flow experiments both inertial and viscous flow terms are not negligible. The observed experimental breakthrough curves of solute transport have been modeled by the classical onedimensional analytical solution for the advection–dispersion equation (ADE) and the single rate mobile–immobile model (MIM). The former model does not properly fit the first arrival and the tail while the latter, which recognizes the existence of mobile and immobile domains for transport, provides a very decent fit. The carried out experiments show that there exists a pronounced mobile–immobile zone interaction that cannot be neglected and that leads to a non-equilibrium behavior of solute transport. The existence of a non-Darcian flow regime has showed to influence the velocity field in that it gives rise to a delay in solute migration with respect to the predicted value assuming linear flow. Furthermore, the presence of inertial effects enhance non-equilibrium behavior. Instead, the presence of a transitional flow regime seems not to exert influence on the behavior of dispersion. The linear-type relationship found between velocity and dispersion demonstrates that for the range of imposed flow rates and for the elected path the geometrical dispersion dominates the mixing processes along the fracture network.


Introduction
Proper management of groundwater resources requires an understanding of the processes that cause water contamination and affect the remediation of polluted aquifers (Cherubini et al., 2010;Cherubini and Pastore, 2011).
In fractured rock aquifers, open fractures as well as bedding planes or faults give place to preferential flow paths for ground water, contaminants in solution, and free product to reach very quickly an exposure point (Becker and Shapiro, 2003).Recently, research has targeted physical mechanisms that can mitigate fast-path transport by delaying mass en route (Neretnieks, 1980;Haggerty and Gorelick, 1994;Ostensen, 1998;Haggerty et al., 2000;Becker and Shapiro, 2003).

Published by Copernicus Publications on behalf of the European Geosciences Union.
A way for understanding the migration of contaminants in such complex systems is that of analyzing tracer transport, measuring mass-exchange processes in situ, and extracting transport model parameters from the data (Becker and Shapiro, 2003).Testing involves injecting tracers and then monitoring the tracer concentrations as a function of location and/or time (Callahan, 2001).A "breakthrough curve" (BTC) is generated as a function of time at the observation well to estimate the transport parameters.Due to high costs of field transport experiments, in many cases the choice is that of carrying out small-scale (laboratory) experiments to analyze solute transport in representative rock samples.
Moreover laboratory experiments give the advantage of improving the understanding of physical mechanisms under relatively well-controlled conditions, since the dependence of a physical process on different parameters can be tested and modeled.
Obtaining reliable predictions of processes that govern contaminant propagation in fractured aquifers requires modeling the advection and the dispersion of contaminants at the scale of the fracture network.In this context understanding transport in a single fracture is a crucial first step (Boschan et al., 2008).Qian et al. (2011) carried out well-controlled laboratory experiments to investigate flow and transport in a single fracture under non-Darcy flow conditions.Non-Fickian transport was found to dominate with early first arrival and long tails.A mobile-immobile (MIM) model proved to fit both peak and tails of the observed BTCs better than the classical ADE (advection-dispersion equation) model.
On the basis of this experience, in order to describe solute transport under different flow velocities and fracture apertures, Chen et al. (2011) carried out a series of wellcontrolled flow and tracer test experiments on an artificial channeled single fracture (CSF) -a single fracture with contact in certain areas -constructed in the laboratory.The flow condition showed a non-Darcy feature (best described by the nonlinear Forchheimer equation) and the BTCs showed a non-Fickian nature of transport such as early arrival of the peak value, long tailing and multi-peak phenomena.The results of this study showed that the ADE is not adequate to describe the BTCs in a CSF.Sudicky et al. (1985) examined the migration of a nonreactive tracer in layered media under controlled laboratory conditions by conducting multiple tracer tests with a column containing stratified porous media.The experimental results showed strongly dispersed and skewed shape of the breakthrough curves in contrast to the symmetric and weakly dispersed concentration patterns typically associated with homogeneous media.Simulations of the experiments demonstrated that these tailing effects are the result of a transient redistribution of the tracer across the strata by transverse molecular diffusion and that local longitudinal dispersion is only of secondary importance as a spreading process in such systems.These findings were consistent with recent theoretical descriptions of dispersion in stratified aquifers.Starr et al. (1985) carried out reactive tracer tests on the same column and found breakthrough curves that were similar in form to those reported in the preceding study for a nonreactive solute, but were delayed in the time of appearance, had a lower peak concentration, and were more highly dispersed.A mathematical model accounting for longitudinal advection in the sand layer, transverse diffusion in the silt layers, and retardation in both the sand and silt layers gave a good representation of the experimental data, however significant discrepancies existed between the measured and simulated results, with the discrepancies becoming greater at lower velocities.The less satisfactory agreement obtained in the latter study suggested that there is some physical or chemical aspect of the retardation process that was not adequately represented in the model.Callahan et al. (2000) carried out laboratory experiments by means of multiple experimental tracer methods to determine fracture/matrix interactions and dispersion in artificially created fractured rock core of volcanic tuff.They affirmed that matrix diffusion serves to increase the transport time of solutes in dual porosity media by spreading mass away from the advecting region of the fractures.However, they also affirmed that the results of these short-term tests were probably influenced to some degree by smaller-scale processes that should be minimal in field experiments, such as diffusion within the stagnant water in the fractures ("freewater diffusion"), caused by fracture aperture variability, that were more important at small time scales.Because free water diffusion coefficients are larger than matrix diffusion coefficients, this led to an overestimation of the amount of diffusive mass transfer (Callahan et al., 2000).Leven et al. (2005) carried out tracer tests at bench scale on artificially fractured laboratory blocks using port-port connections in such a way as to create matrix-dominated and fracture-dominated ports.Breakthrough curves detected at matrix-dominated port connections were characterized by mainly broad and flat curves in contrast to breakthrough curves recorded at the outlet of direct fracture connections, which showed earlier first arrivals with much sharper and steeper concentration increases.
The authors attributed short breakthrough times to a fast transport of the tracer through the fracture system with a less pronounced interaction with the matrix.The broad and flat breakthrough of tracer was attributed to transport mainly through the matrix with dominant diffusive transport mechanisms.The presence of tails in the BTC curves was attributed to matrix diffusion and to differential advection, i.e. the existence of pathways of varying length through the dimensionality of the flow field (McDermott et al., 2003).Rodrigues et al. (2008) carried out several small-scale (laboratory) tracer tests to analyze advection and dispersion of different solutes in fractured media.
The results obtained with sodium chloride revealed differences between direct and reverse tests, due to its higher density than water.As a consequence, in case of significant openings near the injection hole the solution might sink there initially and be released later by diffusion.In case of no pits near the injection point, the solution could mix with the flowing water more easily and therefore travel faster through the system.
Thus, as the transport was found to depend on the morphology of the areas around the injection holes, the parameters calculated in the tests with sodium chloride did not correspond to the parameters of the flowing water.
On the other hand, the tests with fluorescein or sulforhodamine allowed to obtain parameters characterizing the flow pattern without being affected by exogenous factors because they behaved as inert tracers showing advection-dominated transport with high Péclet numbers.
The present study uses well-controlled laboratory experiments to investigate flow and transport in an artificially fractured laboratory block.Flow in the experiments is nonlinear and is well described by the Forchheimer equation (Cherubini et al., 2012).Non-Fickian transport is found to dominate with early first arrival and long tails.The BTCs of the solute transport are modeled by the conventional ADE, and the single rate (MIM) model.The former poorly describes the behavior of the breakthrough curves while the latter is able to fit the peak value and the tail.

Flow models
Generally the model used to describe fluid flow in fractured media is the local cubic law, which adapts Darcy's law under the assumption of ideal fractures with flat, smooth and parallel walls with infinite lengths, together with laminar flow, incompressible fluid and confined configuration.Different studies in literature show that in real rock fractures a nonlinear flow behavior is easy to occur (Kolditz, 2001;Javadi et al., 2010).A flow model commonly used to represent non-Darcy flow behavior is the Forchheimer law, which includes a quadratic term of velocity to represent the inertial effect: where h (L) is the hydraulic head, x (L) is the spatial coordinate along the direction of the flow, v (LT −1 ) is the flow velocity, a (TL −1 ) and b (T 2 L −2 ) are the linear and inertial coefficients respectively.A general Darcian-like relationship can be used to describe nonlinear flow regimes: K eff (LT −1 ) represents the effective hydraulic conductivity function of hydraulic gradient.According to Forchheimer's law, the effective hydraulic conductivity can be written as . (3)

Solute transport models
One of the most widely used solute transport models in field applications reported in literature is the ADE.The ADE model is based on Fick's law, which assumes that the dispersive mass flux is proportional to the first-order spatial derivative of concentration (Bear, 1972).The mathematical formulation of the ADE model for nonreactive solute transport can be expressed as follows: where t (T) is the time, x (L) is the spatial coordinate along the direction of the flow, c (ML −3 ) is the solute concentration, v (LT −1 ) is the average flow velocity and D (L 2 T −1 ) is the dispersion coefficient.In complex, fractured media the latter depends mainly on two processes: Aris-Taylor dispersion (D T ) (Aris, 1956) due to the combined action of convection and radial molecular diffusion (Dullien, 1992) and geometrical dispersion (D G ) due to the roughness and/or aperture variation of fractures (Boschan et al., 2008).For a fracture characterized by two flat parallel walls geometrical dispersion should be equal to zero and dispersion processes are represented only by Aris-Taylor dispersion (Auriault, 1995) by the following expression: where w (L) is fracture aperture, D m (L 2 T −1 ) is the molecular diffusion.Theoretically the above equation is valid only beyond a critical travel time, which corresponds to the minimum duration needed for a particle to experience the whole cross-sectional parabolic profile of velocities across the fracture aperture (Bodin et al., 2007).The critical travel time is proportional to a characteristic time τ c of transverse diffusion equal to In other words the solute must travel over a minimum distance before the Aris-Taylor dispersion regime is completely established.For t < τ c the dispersion coefficient is time dependent (Berkowitz and Zhou, 1996).Geometrical dispersion is nonzero only for fractures with rough walls and/or with variable apertures and reflects the influence of spatial variation of flow velocity in the plane of fractures.The geometrical dispersion resulting from heterogeneity along the fracture plane varies linearly with the mean velocity: where α L (L) is the geometrical dispersivity coefficient.The Péclet number is defined as It can be used to distinguish different regimes in variable aperture fractures: molecular diffusion, geometric dispersion and Aris-Taylor dispersion.The geometric dispersion regime corresponds to the range of P e where velocity variations in the plane of fractures dominate the mixing process.This means that for low values of α L the dispersion passes directly from the molecular diffusion regime to the Aris-Taylor dispersion regime whereas for high values of α L there exists a large range of Pe in which geometrical dispersion dominates.The Péclet number is also defined as Where L (L) is the characteristic length of the domain.It represents the relative effect of advective compared to dispersive solute transport.At high Péclet numbers P e 1 advection dominates solute transport processes; while at low Péclet numbers P e < 1 dispersion/diffusion dominates.
In a typical tracer injection experiment a mass M 0 (M) of the tracer is injected instantaneously at time zero at the origin of the domain (x = 0).The initial condition is given by Where ω (L 2 ) is the cross-sectional area, δ(x) (L −1 ) is the Dirac delta, which is equal to 1 when x is equal to zero and is 0 otherwise.In addition, it is assumed that a first-type boundary condition exists at the outflow boundary: The solution of Eq. ( 4) for the specified initial and boundary conditions is given by (Crank, 1956) Generally, heterogeneous media such as fractured media present non-ideal transport of solutes that lead to asymmetrical BTCs.Dual porosity models or the MIM take into account non-ideal transport behavior caused by the presence of immobile zones.MIM model assumes that the net mass transfer from the main flow field (mobile domain) to the stagnation zones (immobile domain) is proportional to the difference in concentration between the mobile and immobile domains.The mathematical formulation of the MIM for nonreactive solute transport is usually given as follows: where c m and c im are the cross-sectional averaged solute concentrations respectively in the mobile and immobile domain, α (T −1 ) is the mass exchange coefficient, and β [-] is the mobile water fraction.For a nonreactive solute β is equivalent to the ratio between the immobile and mobile cross-sectional area (-).
The Damköhler (Da) number can be used to evaluate the behavior of the MIM model, but it has been showed (Wagner andHarvey, 1997, 2001) that there are limitations to a tracer test's ability to estimate the exchange parameters αand β.
This dimensionless number can be expressed as (Wagner and Harvey, 2001) At high values of Da the solute concentrations in mobile and immobile domains are in equilibrium and MIM tends to the ADE model.In this case the effect of the exchange is difficult to identify.At very low values of Da the mass transfer is absent or very slow and a dual domain performs as a single domain.Only a small amount of tracer interacts with the immobile zones, thus the exchange effect is small and difficult to identify.In between these values of Da the mass transfer is controlled by a first order kinetic process depending on the concentration gradient between the mobile and the immobile domains.Bahr and Rubin (1987) demonstrated the use of this dimensionless number for identifying those cases where nonequilibrium transport cannot be distinguished from equilibrium transport.Da can be viewed as the ratio of the transport time (L/v) to the exchange time equal to the reciprocal of the exchange coefficient (1/α).On the basis of these considerations it can be stated that the effect of the immobile zones is stronger when the exchange time is of the same order of magnitude as the transport time.
In analogous manner for the ADE model the solution of system Eq.( 13) describing one-dimensional nonreactive solute transport in an infinite domain for instantaneous pulse of solute injected at time zero at the origin is given by (Goltz and Roberts, 1986) with where I 1 represents the modified Bessel function of order 1.
In order to simplify the calculation of the analytical solutions, a unit length for the one-dimensional domain has been assumed.It is therefore necessary to normalize the parameters v and D so they both have units of T −1 .To obtain the normalized values of v and D, it is simply necessary to divide them by L and L 2 respectively.

Convolution solutions for variable boundary conditions
The transport solution presented in the previous sections allows the determination and the prediction of BTCs at a specified distance from the inlet boundary.This solution assumes instantaneous pulse injection input condition.In many experimental systems the input boundary condition could be different or the transport occurs through regions with distinctly different properties (Berkowitz et al., 2001).Generally, at field scale both the injection time and the residence time of the solute within the probe are negligible compared to the residence time of solute in the aquifer.Instead, at laboratory scale these assumptions could not be valid because the residence times of the solute in the medium and in the probe are of the same order of magnitude.Convolution techniques can be used to overcome these problems.
Given the BTC curves for the medium C(t) and for the probe S(t) corresponding to the analytical solutions for a pulse input, the convolution of C(t) and S(t) is formally defined as where W (t) is the resulting breakthrough curve recorded by the probe.

Experimental setup
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 in southeastern Italy (Cherubini et al., 2012).The fracture network has been made artificially through 5 kg mallet blows (Fig. 1a).The characteristics of each fracture have been determined from its trace on block surface.The fissured system and the fracture aperture on the block surfaces have been recorded with a high resolution digital camera.Subsequently the images have been scaled and rectified using Perspective Rectifier software (www.rectifiersoft.com)and calibrated on the basis of manual measurements carried out by means of a caliper.Profiles of discontinuities and aperture measurements have been extracted from the recorded images using the 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.
The surface of the block sample has been sealed with transparent epoxy resin (Leven et al., 2004) (Fig. 1b).A hole of 1 cm diameter has been opened for each discontinuity along the boundary of the block by means of a percussion drill (Fig. 1c).Inside each hole an hexagonal bushing of 1/4 M-3/8 F has been placed and welded to the block by means of rapid-hardening epoxy resin (Fig. 1d).
In the present paper the study of flow and transport dynamics regards only a single path.Figure 2a shows the mechanical aperture distribution obtained from 13688 measurements and Fig. 2b shows the three-dimensional geometry of the selected path.The average cross-sectional area of the path is equal to ω = 0.9932 whereas the average path length is equal to L = 0.7531 m.
The sealed block sample is connected to a hydraulic circuit (Fig. 3).Water inside the block flows according to the hydraulic head difference between the upstream tank connected to the inlet port and the downstream tank connected to the outlet port.The upstream and downstream tanks have the same dimensions and are of cylindrical shape with a circular cross section.The instantaneous flow rate that flows across the block is measured by an ultrasonic velocimeter (DOP3000 by Signal Processing).Beside the inlet port a syringe for instantaneous injection of a conservative tracer (sodium chloride, NaCl) is placed, while at the outlet port there is a flow cell in which a probe can be positioned.The probe is a multi-parametric instrument (Idronaut Ocean Seven 304 CTD logger) with frequency of sampling of 8 Hz for instantaneous measurement of pressure (dbar), temperature ( • C) and electric conductivity (µS cm −1 ) respectively with resolutions of 0.0015 %, 0.0006 • C and 0.1 µS cm −1 .

Flow tests
The analysis of flow dynamics through the selected path regards the observation of water flow from the upstream tank to the flow cell with a circular cross section of 0.1963 and 1.28 × 10 −4 m 2 respectively.
Initially at time t 0 , the valves a and b are closed and the hydrostatic head in the flow cell is equal to h 0 .The experiment begins with the opening of 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 experiment procedure is repeated changing the hydraulic head of the upstream tank h c .The time t (t 1 − t 0 ) required to fill the flow cell from h 0 to h 1 has been registered.
Given that the capacity of the upstream tank is much higher than that of the flow cell it is reasonable to assume that during the experiments the level of the upstream tank (h c ) remains constant.Under this hypothesis the flow inside the system is governed by the equation: where S 1 (L 2 ) and h (L) are respectively the section area and the hydraulic head of the flow cell; h c (L) is the hydraulic head of upstream tank, and ( h) represents the hydraulic conductance term representative of both hydraulic circuit and the selected path.Hydraulic loss within the single hydraulic circuit can be expressed according to Chézy's law as where R c (T 2 L −5 ) is a characteristic coefficient related to the roughness, section and length of the tubes of the hy-draulic circuit, h (L) is the hydraulic head difference, and Q (L 3 T −1 ) is the flow rate.Whereas, only for the sealed block, the h-Q relationship can be expressed by means of a discrete form of Forchheimer's law: where A (TL −2 ) and B (T 2 L −5 ) are the linear and nonlinear hydraulic loss coefficients respectively and are related to the roughness, aperture, lengths and shape of the selected path in the fractured medium.
Combining Eqs. ( 19) and ( 20) the hydraulic conductance term of the whole hydraulic system assumes the following expression: The average flow rate Q can be estimated by means of the volumetric method: Whereas the average hydraulic head difference h is given by In correspondence of the average flow rate and head difference is it possible to evaluate the average hydraulic conductance as The inverse of ( h) represents the average resistance to flow.Substituting Eq. ( 21) into Eq.( 18) and integrating the latter from t = t 0 to t = t 1 with the initial condition h = h 0 the following expression is obtained for Forchheimer's law: By fitting the experimental relation between t = t 1 −t 0 and h c it is possible to obtain an estimation of the coefficients A and B.

Tracer tests
The study of solute transport dynamics through the selected path has been carried out by means of a tracer test using sodium chloride.Initially a hydraulic head difference between the upstream tank and downstream tank is imposed.At t = 0 the valve a is closed and the hydrostatic head inside the block is equal to the downstream tank.At t = 10 s the valve a is opened while at time t = 60 s a mass of solute equal to 5 × 10 −4 kg is injected into the inlet port through a syringe.The source release time (1 s) is very small therefore the instantaneous source assumption can be considered valid.
In correspondence of the flow cell in which the multiparametric probe is located it is possible to measure the tracer breakthrough curve and the hydraulic head; in the meanwhile the flow rate entering the system is measured by means of an ultrasonic velocimeter.For different flow rates a BTC curve can be recorded at the outlet port.
In order to determine the two parameters of the ADE model (v and D) and the four parameters of the MIM model (v, D, α and β) it is crucial to estimate the function S(t) that appears in Eq. ( 17).The latter has been evaluated through different tests carried out on the probe subject to pulse injection varying the flow rate.Once known S(t) it is possible to obtain an estimation of ADE and MIM model parameters by a fitting procedure between experimental data and theoretical BTC curves evaluated through the analytical solutions.

Calibration of experimental apparatus
As regards the analysis of flow behavior, in order to estimate R c several tests have been conducted only on the hydraulic circuit varying h c in the range of 0.14-0.93m, whereas the average flow rate Q varies in the range   1.77 × 10 −5 -6.80 × 10 −5 m 3 s −1 .The relationship h ct has been fitted by means of Eq. ( 25) with parameters A and B equal to 0 (Fig. 4).The coefficient R c results as R c = 7.10 × 10 −5 s 2 m −5 .
As regards the study of solute transport, the tracer injection device has been connected directly with the flow cell in which the multi-parametric probe is positioned.Several tracer tests have been conducted on this configuration varying the input flow rate in the range of 3.53 × 10 −6 -5.32 × 10 −6 m 3 s −1 .The observed BTCs show an exponential decay function that can be expressed as follows: where Vol (L 3 ) is the volume of the flow cell, C 0 = M 0 /Vol (ML −3 ) is the concentration observed at t = 0. Figure 5 shows the experimental relationship between the observed Qand Q/Vol obtained through the fitting of BTCs.The relationship is linear and the estimated volume of flow cell Vol is 1.237 × 10 −4 m 3 , close to the real volume of flow cell equal to 1.417 × 10 −4 m 3 .

Flow characteristics
Several tests have been conducted for the selected path.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 1.85 × 10 −6 -1.11 × 10 −5 m 3 s −1 .Figure 6 shows the fitting method described in the previous section to estimate the linear and nonlinear terms equal to A = 4.11 × 10 4 (sm −2 ) and B = 6.61 × 10 9 (s 2 m 5 ) respectively.A nonlinear Q-h relationship that is best described by Forchheimer's equation is shown in Fig. 6b.Flow characteristics can be studied through the analysis of linear (AQ) and quadratic terms (BQ 2 ).These terms represent the head losses due to viscous and inertial forces respectively and they can be used to analyze flow conditions.The ratio of nonlinear to linear head losses is the Forchheimer number: Inertial forces dominate over viscous ones at the critical Forchheimer number (F 0 = 1) corresponding in our case to a flow rate equal to Q = 6.21 × 10 −6 m 3 s −1 .
Considering the height of the block equal to 0.08 m the viscous term is larger than the inertial one as long as q < 77×10 −6 m 2 s −1 .This specific discharge is of the same order of magnitude as the one found by Qian et al. (2011) in a filled single fracture (0.250 m height).
Even if inertial forces become dominant at higher flow rates, Fig. 6c and d show a consistent variation of the conductance and of the resistance terms respectively.The average hydraulic conductance in terms of flow rate can be written as   According to the above expression the conductance for Darcian flow can be approximated with the reciprocal of A and is equal to 2.43 × 10 −5 m 2 s −1 .In correspondence of the critical Forchheimer number the conductance is less than 50 %, reaching the value of 1.22 × 10 −5 m 2 s −1 .On the basis of these considerations it could be argued that in the executed flow experiments both inertial and viscous flow terms are not negligible.Generally, in highly fractured aquifers under natural conditions low hydraulic gradients are detectable and inertial characteristics of the flow field can be negligible.Anthropogenic stresses in the aquifer such as pumping and recharge scenarios give rise to high hydraulic gradients and consequently the impact of inertial effects cannot be neglected.Pumping tests with multiple pressure steps should be used for a more accurate identification of the range of validity of the Darcian-type relationship and for the quantification of linear and nonlinear flow relations.

Solute transport
The aim of this paper is to observe if the presence of inertial effect and the variation of flow rate can influence the behavior of solute transport.To do that several tests have been conducted for the selected path in order to observe solute transport patterns varying the flow rate in the range 1.20 × 10 −6 -8.34 × 10 −6 m 3 s −1 .Each obtained experimental BTC curve has been fitted by ADE and MIM models.for 4 × 10 −6 m 3 s −1 < Q < 6 × 10 −6 m 3 s −1 .Figure 8 shows the measured BTCs fitted by ADE and MIM models simultaneously for given average flow rates of 5.98 × 10 −6 , 2.54 × 10 −6 , 1.84 × 10 −6 and 1.32 × 10 −6 m 3 s −1 .These results demonstrate that the MIM model fits adequately the observed BTCs.However, especially for low values of flow rate the MIM model does not fit precisely the experimental BTCs, in fact a "late peak" is observed, probably due to the presence of a second "delayed" flowpath.
In Fig. 9 the ratio of linear and nonlinear losses to total head loss has been plotted as function of flow rate.On the same diagram the experimental relationships of flow rate and normalized velocity for ADE and MIM models have been superimposed with ADE underestimating the velocity with respect to MIM.In the transition between viscous and inertial regime a change in the slope can be evidenced, which means the setting up of a different behavior.In physical terms this means that the diagram of velocity profile is flattened because of inertial forces prevailing on viscous ones.Therefore the presence of the transitional flow regime leads to a delay on solute transport with respect to the values that can be obtained under the assumption of a linear flow field.
Figure 10  time vs. flow rate can be evidenced and as a consequence the non-equilibrium behavior becomes stronger.
In Fig. 11 the first-order mass exchange term α is plotted as a function of the normalized velocity (v/L).For low normalized velocity values α is constant, while in correspondence of a critical value it starts increasing with a potential function.For the experiment conducted the masstransfer coefficient is dependent on flow velocity.This behavior is attributable to higher flow velocities that enhance the exchange between the main flow paths and the diffusion dominated zone.Several authors have observed variation of the mass-transfer coefficient between mobile and immobile water regions with pore-water velocity (van Genuchten and Wierenga, 1977;Nkedi-Kizza et al., 1984;De Smedt and Wierenga, 1984;De Smedt et al., 1986;Schulin et al., 1987).The increase in α with increasing pore-water velocity is attributed to higher mixing in the mobile phase at high porewater velocities (De Smedt and Wierenga, 1984) or to shorter diffusion path lengths as a result of a decrease in the amount of immobile water (van Genuchten and Wierenga, 1977).
The mobile water fraction β coefficient exhibits a constant mean value of 0.56 and a weak trend growth as velocity increases.This means that the mobile and immobile domains present the same size and confirm the presence of non-equilibrium solute transport behavior.
As regards the Damnkhöler number Da estimated by Eq. ( 14) a linear trend as velocity increases has been detected, its values are in the range 0.39-1.61.According to Wagner and Harvey (1997) this means that the exchange transport parameters found for the MIM model are reliable because Da is in the order of magnitude 0.1-1.0.This demonstrates that there exists a pronounced mobileimmobile zone interaction that cannot be neglected and that leads to a non-equilibrium behavior of solute transport.
Moreover, while the presence of a transitional flow regime affects the velocity field, it seems not to exert influence on the behavior of dispersion.In Fig. 12 a linear relationship is shown between velocity and dispersion both for ADE and MIM models, with ADE overestimating the dispersion with respect to MIM, as found by Qian et al. (2011).These experimental results suggest that geometrical dispersion dominates the effects of Aris-Taylor dispersion.This can be attributable to two possible mixing processes: the presence of a transient regime that precedes stationary dispersion and the presence of complex paths that enhance the geometrical dispersion and do not allow the development of Aris-Taylor dispersion.
In order to get the Aris-Taylor dispersion regime completely established the solute must travel over a minimum distance much larger than the product of velocity and the characteristic time of transverse diffusion (Eq.6).In our case considering a value for molecular diffusion of D m ≈ 1.85 × 10 −9 m 2 s −1 , the ratio between the minimum distance and the length of active path is in the range 2.88-19.61,therefore according to Berkowitz and Zou (1996) the Aris-Taylor dispersion coefficient is time dependent.In the case study, due to the complex, tortuous topology of the fracture network there may not be an opportunity for Aris-Taylor dispersion to develop.The obtained results need to be compared with previous studies carried out in limestone, granite and welded tuff fractures (Kumar et al., 1995;Yeo et al., 1998;Wan et al., 2000) that found the dominance of Aris-Taylor dispersion.The previously reported experiments (Detwiler et al., 2000) regard a single fracture in which the fraction of contact area between the fracture surfaces was relatively small.The presence of long, simply connected flow paths permitted Taylor dispersion to "develop".Instead in our study a fracture network is present, fracture intersections interrupt the continuity of flow paths between single fractures and give rise to velocity fluctuations that enhance geometrical dispersion.

Conclusions
In this paper controlled laboratory experiments on flow and transport have been carried out at bench scale in an artificially fractured limestone block.For a selected path both hydraulic and pulse tracer tests have been conducted.Hydraulic tests have proved the existence of a nonlinear flow behavior best described by the Forchheimer's law.Transition from viscous dominant regime to inertial dominant regime has been detected.
All experimental BTCs obtained from pulse injection tests varying the flow rate exhibit a pronounced non-equilibrium transport regime especially at high flow rates that cannot be neglected.
The presence of not negligible inertial effects plays an interesting role on the structure of the flow field and can lead to the presence of recirculation zones that enhance the non-equilibrium behavior.These recirculation zones represent low velocity or stagnant regions that represent the physical immobile zones.The fact that in the present study the Damköhler number is close to unity especially for high flow rates and that the exchange coefficient depends on flow velocity could confirm this hypothesis.
Mass-transfer limitation can have a significant impact on the mobility of contaminants in groundwater (Maraqa, 2001) and therefore on the predicted cleanup times.
Berglund and Cvetkovic (1995) proposed an analytical solution for the radial mass flux in a heterogeneous aquifer, which refers to a simple case of aquifer remediation by the pump and treat method using an extraction well.The results showed a large initial drop in concentration followed by a leveling and a very slow gradual decline.The general effects of the rate-limited mass-transfer process are a considerable increase in cleanup times and that cleanup times tend to converge towards constant values, where further increase in pumping rate has no effect in terms of reduced cleanup times.They concluded that the advantages of using high pumping rates are small if the recovery of the contaminants in the immobile phase depends on a mass-transfer process characterized by a large timescale.
Several studies (Maraqa, 2001;Zimmerman, 1998) suggested a correlation between the mass-transfer coefficient and pore-water velocity.The increase in the mass-transfer coefficient with increasing pore-water velocity is attributed to higher mixing in the mobile phase at high pore-water velocities or to shorter diffusion path lengths.
The existence of a non-Darcian flow regime has showed to influence the velocity field in that it gives rise to a delay in solute migration with respect to the predicted value assuming linear flow.
Instead, the presence of a transitional flow regime seems not to exert influence on the behavior of dispersion.The linear-type relationship found between velocity and dispersion demonstrates that for the range of imposed flow rates and for the selected path the geometrical dispersion dominates the mixing process along the fracture network.
Quantification of dispersive transport properties is important when delineating protection and capture zones around pumping wells.The determination of the protection zones is based on travel time analysis for an imaginary groundwater pollutant, assuming advective and dispersive transport.
In the case study, an important feature that merits additional analysis is the fundamental difference between the carried out studies on single fractures (Kumar et al., 1995;Yeo et al., 1998;Wan et al., 2000) in which Taylor dispersion was found to occur and the dominance of geometric dispersion in our complex fracture network.
In most fractures used in previously reported experiments (Detwiler, 2000), the fraction of contact area between the fracture surfaces was relatively small.The resulting long, simply connected flowpaths permitted Taylor dispersion to "develop".In the case study, due to the complex, tortuous topology of the fracture network there may not be an opportunity for Taylor dispersion to develop.Furthermore the presence of tortuosity may display very different geometrical features and can have an impact on different transport mechanisms such as preferential flow paths, retardation zones and limited mass exchange between zones.
Increasing evidence indicates that the advectiondispersion equation has problems in describing solute transport in fractured media.A mobile-immobile model has showed to be a better candidate for a conceptual model (Chen et al., 2010) partially because it has two more fitting parameters than the ADE, and also because it has recognized the physical reality of the existence of both mobile and immobile domains.
In this paper, the MIM model proves to fit adequately the observed BTCs except for low flow rate values where a "late peak" is observed, probably due to the presence of a second "delayed" flowpath.Further developments of this study would deal with investigating the performance of more accurate models (Bodin et al., 2007) to predict tracer transport behavior in a fracture network.Moreover, in order to be able to generalize the relationships between models parameters different tracer tests will be carried out for different paths and for different blocks.
In fractured aquifers, the prolonged BTCs can also be attributed to intermediate storage in zones containing quasiimmobile groundwater, and slow release into active fractures.These phenomena demonstrate that fractured aquifers are not always fast-flushing systems, but contaminants can sometimes remain in immobile fluid regions for long periods (Goldscheider, 2008).The persistence of some pollutants creates a costly, remedial challenge in the subsurface.
Due to the complexity of hydrogeologic conditions in tight formations, setting up a decontamination treatment in fractured rock aquifers is a complex undertaking: traditional remediation processes adopted for porous media cleanup do not work well for the removal of most pollutants.Groundwater extraction and treatment systems are not fully effective in speeding up the aquifer's remediation, because of diffusive processes that tend to retard a contamination plume's advance through a fractured rock aquifer and to substantially increase the difficulty of purging pollutants from it.
Therefore, the traditional remediation techniques have to be combined with new innovative techniques that prove to be effective in speeding up the aquifer's cleanup.
In classic depollution strategies such as Pump & Treat system, an induced-hydraulic gradient much higher than the natural gradient is often necessary to clean up contaminated groundwater, where non-Darcian flow is expected to occur.Under these circumstances, the behavior of solute dispersion needs to be evaluated thoroughly for the full range of flow regimes.
This study motivates further experimental work to increase our understanding of the controlling mechanisms and key parameters for flow and transport in fractured media, and thus make better predictions of the infiltration path of a solute, the location of trapped contaminants, the effect of containment or mobilization techniques for pollutants, and the remediation of contaminated fractured aquifers.
Edited by: G. Fogg

Figure 1
Figure 1.a) artificial discontinuities produced by means of 5 kg mallet blows; b) epoxy resin casting; c) example of a hole along the edges of the discontinuities; d) insertion of hexagonal bushing for the connection to the hydraulic circuit.

Figure 2
Figure 2. a) distribution of mechanical aperture evaluated on the 13688 samples b) three dimensional reconstruction of the fracture network obtained from fracture traces on block.The selected path is highlighted.

Fig. 1 .
Fig. 1.(a) artificial discontinuities produced by means of 5 kg mallet blows; (b) epoxy resin casting; (c) example of a hole along the edges of the discontinuities; (d) insertion of hexagonal bushing for the connection to the hydraulic circuit.

Figure 1 Figure 2
Figure 1.a) artificial discontinuities produced by means of 5 kg mallet blows; b) epoxy resin casting; c) example of a hole 676

Fig. 4 .
Fig. 4. Control head h c vs. time for calibration of hydraulic circuit.Circle represents the experimental values obtained from tests carried out only on the hydraulic circuit.The marked line represents Eq. (22) with A and B equal to 0. 688

Figure 5 .Figure 6 .
Figure 5. Calibration of the probe.Q represents the measured flow rate, Q/Vol is estimated by fitting the BTCs curves o 689

Fig. 5 .
Fig. 5. Calibration of the probe.Q represents the measured flow rate, Q/Vol is estimated by fitting the BTCs curves of the probes.The circle line represents the experimental values, the dashed line represents the linear regression.

Figure 5 .Figure 6 .Fig. 6 .
Figure 5. Calibration of the probe.Q represents the measured flow rate, Q/Vol is estimated by fitting the BTCs curves of the 689

Figure 8 .
Figure 8. Fitting of BTCs for different flow rate values.cQ/M0 represents the normalized flow mass rate (c is the 701

Figure 9 .FigureFig. 9 .
Figure 9. Flow rate vs. the ratio of linear and nonlinear losses to total loss (black line and red line respectively) Fig. 10.Transport time (L/v) (reciprocal of normalized velocity) and exchange time (1/α) (reciprocal of the exchange term) as a function of flow rate.

Fig. 11 .
Fig. 11.Normalized flow velocity vs. first-order mass exchange coefficient α estimated for the MIM model.
lized flow velocity vs. normalized dispersion estimated for ADE and MIM model.