Fluid Flow, Mass, and Heat Transport Laboratory Experiments in Artificially Fractured Rock

The knowledge of flow, mass, and heat transport phenomena in fractured rocks represents an important issue in many situations of science and engineering, such as the groundwater resource management, fractured petroleum reservoirs, nuclear waste disposal, as well as geothermal energy development. Experimental data obtained under controlled conditions such as in laboratory increase the knowledge of the fundamental physics that interest the flow and transport in fracture media and allow us to investigate the behavioral differences between fracture and porous medium. An artificially fractured rock sample of parallelepiped shape (0.60 0.80 0.08 m) has been created, and the flow, mass, and heat transport behavior have been observed. The carried out experiments show the existence of non-Darcian flow regime that cannot be neglected. The latter influences the mass transport behavior giving rise to a delay in mass migration, enhancing the nonequilibrium behavior, whereas the dispersion phenomena seem not be influenced. Heat transport shows a very different behavior compared to mass transport. Convective thermal velocity is lower than mass velocity, whereas thermal dispersion is higher than solute dispersion.


Introduction
In the study of hydrodynamic processes in fractured media, flow is described with a linear relationship between pressure gradient and flow rate [1], valid at low flow regime (Re < 1). For Re > 1, a nonlinear flow can be observed.
In fracture networks, heterogeneity intervenes even in mass and heat transport: due to the variable aperture and heterogeneities of the fracture surfaces, the fluid flow will seek out preferential paths through which mass and heat are transported.
The study of mass and heat transport behavior is based on dual porosity theory. Fractured medium is conceptualized in two domains: connected fractured network where transport processes is governed mainly by advection and rock matrix as well as poorly connected fractures which represent stagnant zones.
The present chapter reports several laboratory experiences at bench scale and their interpretation of flow and transport behavior on artificially fractured rock sample of parallelepiped shape. The experimental results of flow and transport tests in a fractured block at bench scale are interpreted by means of explicit network model (ENM) where the fracture network geometry has been taken into accounts.

Artificial fractured rock sample preparation and characterization
A limestone block of parallelepiped shape (0.6 Â 0.4 Â 0.08 m 3 ) has been recovered from the "Calcare di Altamura" formation located in the Apulian region in Southern Italy. Using a 5-kg mallet blow, a fractured network has been made artificially. The fissured system and their aperture have been recorded with a high-resolution digital camera. The Perspective Rectifier (www.rectifiersoft.com) has been used in order to scale and rectify the picture. The fracture traces and the aperture measurement have been extracted from recorded images through the built-in Scilab Image Processing Toolbox (www.scilab.org) [2] (Figure 1). Once the artificially fractured block sample has been characterized, the surface of the block sample has been sealed with the epoxy resin, and a hole of 1 cm of diameter has been opened on the boundary of the block sample in correspondence of each discontinuity. Finally, an extruded polystyrene panel with thermal conductivity equal to 0.034 Wm À1 K À1 and thickness of 0.05 m has been used to thermally insulate the fractured block sample ( Figure 2).

Experimental setup flow, mass, and heat transport test
The sealed and thermally insulated fractured block sample is connected with a hydraulic circuit. The sketch of the experimental apparatus is shown in Figure 3. Water moves from the upstream tank to the downstream tank and returns to the upstream by means of a pump. An electric boiler with a volume of 10 À2 m 3 has been used to heat the water. The instantaneous flow rate that flows across the block is measured by an ultrasonic velocimeter. Beside the inlet port, a syringe for instantaneous injection of a conservative tracer (NaCl) has been placed, while at the outlet port, there is a flow cell in which a multiparametric instrument can be positioned. Two thermocouples have been placed at the inlet and the outlet of a selected fracture path of the limestone block. A TC-08 thermocouple data logger (Pico Technology) with a sampling rate of 1 s has been connected to the thermocouples.
The study of flow, mass and heat transport regarding the active path highlighted in Figure 3. Figure 4a shows the mechanical aperture distribution obtained from 13,688 measurements and Figure 3b shows the reconstructed three-dimensional geometry of the selected path. The average cross-sectional area of the path is equal to 993 μm whereas the average path length is equal to 0.7531 m [3].

Flow test
The air in the hydraulic circuit and the fractured media has been removed bringing the system to full saturation with the hydrostatic head equal to h 0 . Opening the valve "a," water flows in the system, and the snapshot flow rates are measured by the ultrasonic velocimeters. When the  hydraulic head in the flow cell reaches h 1 , the valve "a" is closed, and opening the valve "b," the hydraulic head in the flow cell is reported to the value h 0 .
The volumetric method can be used in order to estimate the injection flow rate that flows within the fractured block: where S 1 (L) is the cross-sectional area of the flow cell and Δt is the time required to fill the flow cell from h 0 to h 1 . The average flow rate estimated was compared with the snapshot flow rate measured through the ultrasonic velocimeter in order to check the absence of the leaks and losses due to obstruction in the hydraulic system. The experiment is repeated changing the hydraulic head h c (L) of the upstream tank.
Given that the capacity of the upstream tank S 2 (L) is higher than S 1 and the compressibility of the block and the hydraulic circuit is very low, the flow dynamics can be described through the following expressions: where h (L) is the hydraulic head of the downstream flow cell and κ(Δh) is the hydraulic conductance term representative of both the hydraulic circuit and the active fracture network configuration.
For the fractured media, a Δh-Q relationship can be represented through the following polynomial expression: 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 fractures.
The hydraulic head differences Δh are given by The hydraulic conductance term κ(Δh) (L 2 T À1 ) representative of the whole hydraulic system assumes the following expression: The inverse of κ(Δh) represents the average resistance to flow R(Q) (TL À2 ).
Substituting Eq. (5) in Eq. (2) 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 experimental relationship between the time Δt = t 1t 0 , an estimate of parameters A and B can be made.

Mass transport test
Tracer tests by means of sodium chloride have been used in order to study the mass transport dynamics through the selected path. First, a steady-state flow condition has been carried out imposing a hydraulic head difference between the upstream tank and downstream tank. Using a syringe, a mass of sodium chloride of 5 Â 10 À4 kg is injected into the inlet port. The pulse injection assumption can be considered because the source release time is very small (1 s).
The multiparametric probe positioned in the flow cell measures the breakthrough curve (BTC) of the injected sodium chloride and the hydraulic head, in the meanwhile the flow rate is measured using the ultrasonic velocimeter. For different flow rates, a BTC can be registered at the outlet port.

Heat transport test
The study of heat transport dynamics in the selected path takes place by injecting a volume of 10 l of hot water into the block sample. At time t = 0 s, the hydrostatic heads in the block sample and in the downstream tank are equal. At time t = 10 s, valve "a" is opened and water flows into the block sample. At time t = 60 s, the valve "d" is opened, and in the same time the valve "c" is closed. In such a way, the hot water enters in the block sample, and the step temperature function at the inlet port is measured using the first thermocouple. The second thermocouple placed inside the outlet port measures the thermal BTC. Analogous manner to the solute transport of the experiment is repeated varying the flow rate.

Calibration of the experimental apparatus
A calibration procedure is needed in order to eliminate the effect of the experimental apparatus on the measure of the state variables.
As far as the flow experiment, the hydraulic circuit gives rise to a hydraulic loss that must be estimated. Hydraulic losses due to only the hydraulic circuit can be expressed according the Chezy law: where C (L 5/2 T) is the characteristic coefficient related to the roughness, section, and length of the tubes of the hydraulic circuit. Then, C À2 is added to the parameter B in Eq. 6.
The parameter C is estimated conducting several hydraulic tests only on the hydraulic circuit in the range 0.14-0.93 m, corresponding to the average flow rate in the range of 1.77 Â 10 À5 -6.80 Â 10 À5 m 3 s À1 . The relationship h c -Δt has been fitted using Eq. (6) with the parameters A and B equal to 0. The coefficient C result equals to 0.14 Â 10 5 m 5/2 s À2 .
In order to overcome the problem linked with the fact that the residence time of the solute in the probe and in the fracture network are the same order of magnitude, convolution technique has been used.
The BTC recorded by the probe w(t) is represented by the convolution product between the BTC of the fracture network c(t) and the BTC of the probe s(t): The tracer injection device has been directly with the flow cell in which the multiparametric prove is positioned. Several tracer tests have been conducted on this configuration varying the input flow rate in the range 3.53 Â 10 À6 m3 s À1 to 5.32 Â 10 À6 m 3 s À1 . The observed BTCs show an exponential decay function like where Vol (L 3 ) is the volume of flow cell and c 0 = M 0 /Vol (ML À3 ) is the concentration observed at t = 0. The observed BTCs have been fitted using Eq. (9), and the volume of the flow cell Vol has been estimated resulting equally to 1.237 Â 10 À4 m 3 closing the real volume of flow cell equal to 1.417 Â 10 À4 m 3 .
For the heat transport tests, no correction is required; thermocouples at the inlet port and outlet port measure directly the thermal BTCs of the media.

Explicit network model
The 2D explicit network model considers the fracture network geometry permitting a more accurate estimation of flow and transport dynamics. The single fractures are represented as one-dimensional pipe elements forming a 2D pipe network [4].
Assuming that a single fracture (SF) j can be represented by a one-dimensional pipe element, the relationship between the head loss Δh j (L) and the flow rate Q j (L 3 T À1 ) can be written in finite terms on the basis of Forchheimer model [5]: where l j (L) represents the length of fracture and a (TL À3 ) and b (T 2 L À6 ) represent the Forchheimer parameters in finite terms. The term in the square brackets represents the resistance to flow R j (Q j ) (TL À2 ) of fracture j.
The flow field can be obtained in a simple way using Kirchhoff's laws. The flow rate Q j crossing the generic fracture j along a parallel circuit is as follows: where ΣQ (LT À3 ) is the sum of the discharge flow evaluated for the fracture intersection located at the inlet bond of j fracture.
Regarding the mass and transport phenomena, once known as the flow field in the fracture network, the probability density function (PDF) of the residence time at the generic node can be obtained as the sum of the PDFs of each elementary path that reach the node. The latter can be expressed as the convolution product of PDF of each individual fracture belonging to the elementary path. Using the convolution theorem, the PDF at generic node Γ (t) in the network can be represented as follows: where L is the Laplace transform operator, N ep is the number of elementary paths, n f,j is the number of the fractures in i th elementary path, Γ j s ð Þ represents the PDF of the j th SF in the Laplace space, s is the integral variable of the Laplace transform, and P j is the probability of the particle transition at the inlet bond of each individual SF. The rules for particle transition through the fracture intersections play an important role in mass and heat transport. The simplest rule is represented by the "perfect mixing model" in which the mass sharing is proportional to the relative discharge flow rates. Based on the perfect mixing model, the probability of the particle transition crossing the SF is The term in the square brackets in Eq. (4) corresponds to P j .
The PDF represents the density distribution of the time that the mass or heat spends in the fracture network.
In order to know Γ j s ð Þ, a transport model can be assumed and consequently the transport parameters of each SF.
Ref. [6] presented an analytical solution for solute transport in a semi-infinite single fracture embedded in a porous rock matrix. Given that the governing equation of heat and mass transport highlights similarities between the two processes, the analytical solution for solute transport can be used also for heat transport.
According of Tang's solution, the PDF in Laplace space can be expressed as In terms of mass transport, the coefficients υ, A, and β 2 are expressed as follows: where u f (LT À1 ) is the convective velocity, D f (L 2 T À1 ) is the dispersion, θ m is the matrix porosity, D e (L 2 T À1 ) represents the effective molecular diffusion, and δ (L) is the thickness of the boundary layer. For small fractures, δ became the aperture w f (L) of the SF.
In terms of heat transport, the coefficient of the PDF is where θ ¼ r m c m =r w c w and D e ¼ k e = r w c w ð Þ. r m (ML À3 ), c m (ML À3 T À3 K À1 ) and r w (ML À3 ), c w (ML À3 T À3 K À1 ) represent the density and the specific heat capacity of the matrix and water, respectively.
Definitely, the BTC describing for the mass and heat transport in the fracture network as function of time at the generic node, using the convolution theorem, can be obtained as follows: where c 0 (ML À3 ) is the initial concentration, c inj (ML À3 ) is the injection function, and * is the convolution operator. Abate et al. (2006) algorithm is used to perform the inverse Laplace transform.
For heat transport dynamics, the temperature BTC at the generic node of the fracture network is where T 0 (K) is the initial temperature and T inj (K) is the temperature injection function. Three characteristic time scales can be defined as where L (L) represents the characteristic length and t u (T), t d (T), and t e (T) are the characteristic time scales of convective transport, dispersive transport, and loss of the mass or heat into the surrounding matrix [7].
Peclet number (Pe) can be defined as Transport processes are convectively dominated when Peclet number is high; on the contrary, when Peclet number is low, they are dominated by dispersion.
In order to investigate the influence of diffusion in the matrix on convective phenomena, the Damkohler number can be used. Da can be defined as where α (T À1 ) is the exchange rate coefficient corresponding to The exchange rate coefficient α is equivalent to the inverse of t e . When t e reaches t u (Da = 1), concentration and temperature distribution profiles are subject to dual porosity effect characterized by a long tail.
When t e ≫ t u (Da ≪ 1Þ, the fracture-matrix exchange is very slow, and it does not influence mass or heat propagation. On the contrary, when t e ≪ t u (Da ≫ 1Þ, the fracture-matrix exchange is rapid, there is an instantaneous equilibrium between fracture and matrix, and they have the same concentration or temperature. These two circumstances close the standard advective-dispersive transport equation.

Flow
Different flow tests that have been carried out vary the control head h c in the range of 0.17-1.37 m observing the average flow rates in the range of 1.85e-6-1.11e-5 m 3 /s. The linear and nonlinear terms are equal to A = 4.11e4 and B = 6.61e9, respectively. In Figure 5 the fitting results are reported.
Once the linear and nonlinear terms are known, the value of the resistance to flow can be determined as function of the injection flow rate.
The flow field in each single fracture can be estimated in analytical way through the application of Kirchhoff's laws. In Figure 6 a sketch of the 2D pipe conceptualization of the fracture network is reported.
Constant Forchheimer parameter for the whole fracture network has been assumed. They can be derived matching the experimental resistance to flow evaluated experimentally as the inverse of the conductance Eq. (5) with the resistance to flow evaluated as with Q 1 evaluated using the following iterative equation and Q 2 evaluated as The Forchheimer parameters have been estimated, and they are, respectively, equal to a = 7.345 Â 10 4 sm À3 and b = 11.65 Â 10 9 s 2 m À6 .
In order to analyze the experimental results, two dimensionless numbers can be evaluated: the Reynolds number and the Forchheimer number.
Reynolds number (Re) 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 (L) is the hydraulic diameter. For the fracture having small aperture with respect to its height, the hydraulic diameter is equal to D h = w f /2. Re can be reformulated as where d is the thickness of the fractured block.
Forchheimer number (F o ) represents the ratio of nonlinear to linear hydraulic gradient contribution: Figure 6. Two-dimensional pipe network conceptualization of the fracture network of the fractured rock block in Figure 1. Q 0 is the injection flow rate, and Q 1 and Q 2 are the flow rates that are flowing in the parallel branches 6 and 3-4-5, respectively (from Cherubini et al. [9]).
The Forchheimer number can be used for evaluate non-Darcian flow. Inertial effects dominate over viscous effects at the critical Forchheimer number (F o > 1) [8].
Reynolds number 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. High microscopic Reynolds number does not necessarily imply non-Darcian flow. Forchheimer number takes into account both velocity and structure of the medium because the nonlinear term is structure dependent. The linear 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 the nonlinear term will be small, and the Forchheimer number will remain small until the Reynolds number is large. Instead, both the nonlinear term and the Forchheimer number will be large if the structure of the medium is such that microscopic inertial effects can be expected.
For the range of injected flow rate, investigated Re is in the range 11.56-69.37, and Fo are in the range 0.29-1.76. Inertial forces dominate the viscous ones when Fo = 1 corresponding to a flow rate equals to Q crit = 6.30 Â 10 À6 m 3 s À1 . Figure 7. Velocity u f (mÁs À1 ) as a function of the injection flow rate Q 0 (m 3 s À1 ) for ENM with Tang's solution for both mass transport and heat transport (from Cherubini [9]).
The term in square bracket in Eq. (26) is equal to the probability of the particle transition evaluated for the branch 6 as function of Q 0 . When the flow injection increases, the probability of the particle transition decreases. This means that when the injection flow rate increases, the resistance to flow of branch 6 increases faster than the resistance to flow of the branches 3-5.

Transport
The behavior of mass and heat transport has been compared varying the injection flow rates.
The observed heat and mass of BTCs for different flow rates have been individually fitted using the ENM approach. For simplicity the transport parameters u f , D f , and α are assumed equal for all branches of the fracture network.
The experimental BTCs are fitted using Eq. (19) and Eq. (20) for mass and heat transport, respectively.
The determination coefficient (r 2 ) and the root-mean-square error (RMSE) have been used in order to evaluate the goodness of fit. Figure 8. Dispersion D f (mÁs À2 ) as a function of velocity u f (mÁs À1 ) for ENM with Tang's solution for both mass transport and heat transport (from Cherubini [9]).
The results highlight that the estimated convective velocities u f for heat transport are lower than for mass transport, whereas the estimated dispersion D f for heat transport is higher than for mass transport.
A different behavior has been observed for the transfer rate coefficient for mass and heat transport. For mass transport, it can be neglected relatively to the convective velocity, whereas for heat transport, the transfer rate coefficient reaches the convective velocity, and having a characteristic length equal to L = 0.601 m, a dual porosity effect is evident.
In Figure 7 the relationship between u f and Q 0 is reported. Regarding mass transport, experiments for values of Q 0 higher than 4 Â 10 À6 m 3 s À1 u f increase less rapidly. This behavior was due to the presence of inertial forces that gave rise to a retardation effect on solute transport.
A very different behavior is observed for heat transport. Heat convective velocity does not seem to be influenced by the presence of the inertial force, whereas u f is influenced by fracturematrix exchange phenomena resulting in a significant retardation effect. Figure 8 reported the dispersion coefficient D f . Linear relationship between u f and D f suggests that inertial forces did not exert any effect on dispersion and that geometrical dispersion dominates the Aris-Taylor dispersion. Even if heat convective velocity is lower than solute advective velocity, the longitudinal thermal dispersivity assumes higher values than the longitudinal solute dispersivity. Figure 9. Transfer coefficient α (s À1 ) as function of velocity u f (mÁs À1 ) for both mass transport and heat transport (from Cherubini et al. [9]). Figure 9 shows the exchange rate coefficient α as function of the convective velocity u f for both mass transport and heat transport. Regarding the mass transport, the fracture-matrix interaction does not exert any influence. The observed non-Fickian behavior is attributable to the presence of the inertial force and secondary path. Figure 10 shows for mass and heat transport the relationship between Pe and Q 0 . For mass transport Pe increases as Q 0 increases, reaching a constant values Pe = 7.5. For heat transport Pe has a constant trend lower than the unit. Figure 11 shows for mass and heat transport the relationship between Da and Q 0 . For mass transport Da presents a very low value. The Fickian model can be assumed to represent mass transport in each single fracture.
For heat transport Da reaches the unit. It presents a downward trend as Q 0 increases switching from higher to lower values than unit.
Even if Da presents a downward trend as Q 0 increases, when the latter exceeds Q crit , a weak growth trend for Da is detected that however assumes values lower than the unit. Figure 10. Peclet number as function of injection flow rate Q 0 (m 3 s À1 ) for both mass transport and heat transport (from Cherubini et al. [9]).

Conclusion
Laboratory experiments on the observation of flow, mass, and heat transport in a fractured rock sample have been carried out. The parameters that control flow, mass, and heat transport have been estimated using the ENM model. The explicit network model is a useful technique to describe flow, mass, and heat dynamics in fractured media, reducing 2D and/or 3D fracture geometry to a network of 1D pipe elements. However, the real case study can be impracticable because the ENM model approach needs the full knowledge of the fracture network. In order to overcome these difficulties, the ENM can be coupled with equivalent porous media model to represent the bigger fracture which represents the main pathways for flow, mass, and heat dynamics.
Regarding the flow process, the experiments highlighted the dependency between the hydraulic conductivity and the specific discharge in fractured media. Lous and Maini (1970) and Elsworth and Doe (1986) found an underestimation errors of several orders of magnitude in the estimation of transmissivity using non-Darcian constant head data derived by packer tests Figure 11. Da number as function of injection flow rate Q 0 (m 3 s À1 ) for both mass transport and heat transport (from Cherubini et al. [9]).
in fractured rock. The experiments show an underestimated error of 46.59% (average value) of the Darcian flow hydraulic transmissivity.
The flow experiments demonstrate that the flow behavior in the fracture network can be described using the Forchheimer law. Due to the presence of nonlinear flow regime, the probabilities of water distribution between the main and secondary paths are the function of the injection flow rate decreasing as injection flow rate increases.
Several dissimilarities have been detected between heat transport and mass transport. The main discrepancies regarding the transport parameters are as follows: convective solute velocity is higher than heat velocity, whereas solute dispersion is higher than heat dispersion, and the solute exchange rate is negligible compared with the solute velocity, whereas for heat transport, the exchange rate is comparable with the heat velocity giving rise to a very strong dual-porosity effect which delays the heat propagation. On the contrary for the solute transport, the non-Fickian behavior observed in solute BTCs is attributable to the presence of the secondary path combined with the nonlinear flow regime.