Interactive comment on “Laboratory experimental investigation of heat transport in fractured media” by Claudia Cherubini et al

Abstract. Low enthalpy geothermal energy is a renewable resource that is still underexploited nowadays in relation to its potential for development in society worldwide. Most of its applications have already been investigated, such as heating and cooling of private and public buildings, road defrosting, cooling of industrial processes, food drying systems or desalination. Geothermal power development is a long, risky and expensive process. It basically consists of successive development stages aimed at locating the resources (exploration), confirming the power generating capacity of the reservoir (confirmation) and building the power plant and associated structures (site development). Different factors intervene in influencing the length, difficulty and materials required for these phases, thereby affecting their cost. One of the major limitations related to the installation of low enthalpy geothermal power plants regards the initial development steps that are risky and the upfront capital costs that are huge. Most of the total cost of geothermal power is related to the reimbursement of invested capital and associated returns. In order to increase the optimal efficiency of installations which use groundwater as a geothermal resource, flow and heat transport dynamics in aquifers need to be well characterized. Especially in fractured rock aquifers these processes represent critical elements that are not well known. Therefore there is a tendency to oversize geothermal plants. In the literature there are very few studies on heat transport, especially on fractured media. This study is aimed at deepening the understanding of this topic through heat transport experiments in fractured networks and their interpretation. Heat transfer tests have been carried out on the experimental apparatus previously employed to perform flow and tracer transport experiments, which has been modified in order to analyze heat transport dynamics in a network of fractures. In order to model the obtained thermal breakthrough curves, the Explicit Network Model (ENM) has been used, which is based on an adaptation of Tang's solution for the transport of the solutes in a semi-infinite single fracture embedded in a porous matrix. Parameter estimation, time moment analysis, tailing character and other dimensionless parameters have permitted a better understanding of the dynamics of heat transport and the efficiency of heat exchange between the fractures and the matrix. The results have been compared with the previous experimental studies on solute transport.

Abstract.Low enthalpy geothermal energy is a renewable resource that is still underexploited nowadays in relation to its potential for development in society worldwide.Most of its applications have already been investigated, such as heating and cooling of private and public buildings, road defrosting, cooling of industrial processes, food drying systems or desalination.
Geothermal power development is a long, risky and expensive process.It basically consists of successive development stages aimed at locating the resources (exploration), confirming the power generating capacity of the reservoir (confirmation) and building the power plant and associated structures (site development).Different factors intervene in influencing the length, difficulty and materials required for these phases, thereby affecting their cost.
One of the major limitations related to the installation of low enthalpy geothermal power plants regards the initial development steps that are risky and the upfront capital costs that are huge.
Most of the total cost of geothermal power is related to the reimbursement of invested capital and associated returns.
In order to increase the optimal efficiency of installations which use groundwater as a geothermal resource, flow and heat transport dynamics in aquifers need to be well characterized.Especially in fractured rock aquifers these processes represent critical elements that are not well known.Therefore there is a tendency to oversize geothermal plants.
In the literature there are very few studies on heat transport, especially on fractured media.This study is aimed at deepening the understanding of this topic through heat transport experiments in fractured networks and their interpretation.
Heat transfer tests have been carried out on the experimental apparatus previously employed to perform flow and tracer transport experiments, which has been modified in order to analyze heat transport dynamics in a network of fractures.In order to model the obtained thermal breakthrough curves, the Explicit Network Model (ENM) has been used, which is based on an adaptation of Tang's solution for the transport of the solutes in a semi-infinite single fracture embedded in a porous matrix.
Parameter estimation, time moment analysis, tailing character and other dimensionless parameters have permitted a better understanding of the dynamics of heat transport and the efficiency of heat exchange between the fractures and the matrix.The results have been compared with the previous experimental studies on solute transport.

Introduction
An important role in transport of natural resources or contaminant transport through subsurface systems is given by fractured rocks.Interest in the study of dynamics of heat transport in fractured media has grown in recent years because of the development of a wide range of applications, including geothermal energy harvesting (Gisladottir et al., 2016).
Quantitative geothermal reservoir characterization using tracers is based on different approaches to predicting thermal Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
The characterization and modeling of heat transfer in fractured media are particularly challenging as open and wellconnected fractures can induce highly localized pathways which are orders of magnitude more permeable than the rock matrix (Klepikova et al., 2016;Cherubini and Pastore, 2011).
The study of solute transport in fractured media has recently become a widespread research topic in hydrogeology (Cherubini, 2008;Cherubini et al., 2008Cherubini et al., , 2009Cherubini et al., , 2013d;;Masciopinto et al., 2010), whereas the literature about heat transfer in fractured media is somewhat limited.Hao et al. (2013) developed a dual continuum model for the representation of discrete fractures and the interaction with the surrounding rock matrix in order to give a reliable prediction of the impacts of fracture-matrix interaction on heat transfer in fractured geothermal formations.Moonen et al. (2011) introduced the concept of a cohesive zone which represents a transition zone between the fracture and undamaged material.They proposed a model to adequately represent the influences of fractures or partially damaged material interfaces on heat transfer phenomena.Geiger and Emmanuel (2010) found that matrix permeability plays an important role in thermal retardations and attenuation of thermal signals.At high matrix permeability, poorly connected fractures can contribute to the heat transport, resulting in heterogeneous heat distributions in the whole matrix block.For lower matrix permeability heat transport occurs mainly through fractures that form a fully connected pathway between the inflow and outflow boundaries, which results in highly non-Fourier behavior characterized by early breakthrough and long tailing.
Numerous field observations (Tsang and Neretnieks, 1998) show that flow in fractures is being organized in channels due to the small-scale variations in the fracture aperture.Flow channeling causes dispersion in fractures.Such channels will have a strong influence on the transport characteristics of a fracture, such as, for instance, its thermal exchange area, crucial for geothermal applications (Auradou et al., 2006).Highly channelized flow in fractured geologic systems has been credited with early thermal breakthrough and poor performance of geothermal circulation systems (Hawkins et al., 2012).Lu et al. (2012) conducted experiments of saturated water flow and heat transfer in a regularly fractured granite at meter scale.The experiments indicated that the heat advection due to water flow in vertical fractures nearest to the heat sources played a major role in influencing the spatial distributions and temporal variations of the temperature, impeding heat conduction in the transverse direction; such an effect increased with larger water fluxes in the fractures and decreased with a higher heat source and/or a larger distance of the fracture from the heat source.Neuville et al. (2010) showed that fracture-matrix thermal exchange is highly affected by the fracture wall rough-ness.Natarajan et al. (2010) conducted numerical simulation of thermal transport in a sinusoidal fracture-matrix coupled system.They affirmed that this model presents a different behavior with respect to the classical parallel plate fracturematrix coupled system.The sinusoidal curvature of the fracture provides high thermal diffusion into the rock matrix.
Ouyang (2014) developed a three-equation local thermal non-equilibrium model to predict the effective solid-to-fluid heat transfer coefficient in geothermal system reservoirs.They affirmed that due to the high rock-to-fracture size ratio, the solid thermal resistance effect in the internal rocks cannot be neglected in the effective solid-to-fluid heat transfer coefficient.Furthermore the results of this study show that it is not efficient to extract the thermal energy from the rocks if fracture density is not large enough.
Analytical and semi-analytical approaches have been developed to describe the dynamics of heat transfer in fractured rocks.Such approaches are amenable to the same mathematical treatment as their counterparts developed for mass transport (Martinez et al., 2014).One of these is the analytical solution derived by Tang et al. (1981).
While the equations of solute and thermal transport have the same basic form, the fundamental difference between mass and heat transport is that (1) solutes are transported through the fractures only, whereas heat is transported through both fractures and matrix, and (2) the fracturematrix exchange is large compared with molecular diffusion.This means that the fracture-matrix exchange is more relevant for heat transport than for mass transport.Thus, matrix thermal diffusivity strongly influences the thermal breakthrough curves (BTCs) (Becker and Shapiro, 2003).
Contrarily, since the heat capacity of the solids will retard the advance of the thermal front, the advective transport for heat is slower than for solute transport (Rau et al., 2012).
The quantification of thermal dispersivity in terms of heat transport and its relationship with velocity has not been properly addressed experimentally and has conflicting descriptions in the literature (Ma et al., 2012).
Most studies neglect the hydrodynamic component of thermal dispersion because of thermal diffusion being more efficient than molecular diffusion by several orders of magnitude (Bear, 1972).Analysis of heat transport under natural gradients has commonly neglected hydrodynamic dispersion (e.g., Bredehoeft and Papadopulos, 1965;Domenico and Palciauskas, 1973;Taniguchi et al., 1999;Reiter, 2001;Ferguson et al., 2006).Dispersive heat transport is often assumed to be represented by thermal conductivity and/or to have little influence in models of relatively large systems and modest fluid flow rates (Bear, 1972;Woodbury and Smith, 1985).Some authors suggest that thermal dispersivity enhances the spreading of thermal energy and should therefore be part of the mathematical description of heat transfer in analogy to solute dispersivity (de Marsily, 1986), and have incorporated this term into their models (e.g., Smith and Chapman, 1983;Hopmans et al., 2002;Niswonger and Prudic, 2003).In the same way, other researchers (e.g., Smith and Chapman, 1983;Ronan et al., 1998;Constanz et al., 2002;Su et al., 2004) have included the thermomechanical dispersion tensor representing mechanical mixing caused by unspecified heterogeneities within the porous medium.
By contrast, some other researchers argue that the enhanced thermal spreading is either negligible or can be described simply by increasing the effective diffusivity; thus, the hydrodynamic dispersivity mechanism is inappropriate (Bear, 1972;Bravo et al., 2002;Ingebritsen and Sanford, 1998;Keery et al., 2007).Constantz et al. (2003) and Vandenbohede et al. (2009) found that thermal dispersivity was significantly smaller than the solute dispersivity.Others (de Marsily, 1986;Molina-Giraldo et al., 2011) found that thermal and solute dispersivity was on the same order of magnitude.
Tracer tests of both solute and heat were carried out at Bonnaud, Jura, France (de Marsily, 1986), and the thermal dispersivity and solute dispersivity were found to be of the same order of magnitude.Bear (1972), Ingebritsen andSanford (1998), andHopmans et al. (2002), among others, concluded that the effects of thermal dispersion are negligible compared to conduction and set the former to zero.
However, Hopmans et al. (2002) showed that dispersivity is increasingly important at higher flow water velocities, since it is only then that the thermal dispersion term is of the same order of magnitude or larger than the conductive term.Sauty et al. (1982) suggested that there was a correlation between the apparent thermal conductivity and Darcy velocity; thus, they included the hydrodynamic dispersion term in the advective-conductive modeling.
Other similar formulations of this concept are present in the literature (e.g., Papadopulos and Larson, 1978;Smith and Chapman, 1983;Molson et al., 1992).Such treatments have not explicitly distinguished between macrodispersion, which occurs due to variations in permeability over larger scales, and the components of hydrodynamic dispersion that occur due to variations in velocity at the pore scale.
The present study is aimed at providing a better understanding of heat transfer mechanisms in fractured rocks.Laboratory experiments on mass and heat transport in a fractured rock sample have been carried out in order to analyze the contribution of thermal dispersion in heat propagation processes, the influence of nonlinear flow dynamics on the enhancement of thermal matrix diffusion and finally the optimal conditions for thermal exchange in a fractured network.
Section 1 shows a short review of mass and heat transport in fractured media highlighting what is still unresolved or contrasting in the literature.
In Sect. 2 the theoretical background related to nonlinear flow and solute and heat transport behavior in fractured media has been reported.
A better development of the Explicit Network Model (ENM), based on Tang's solution developed for solute transport in a single semi-infinite fracture inside a porous matrix, has been used for the fitting of the thermal BTCs.The ENM model explicitly takes the fracture network geometry into account and therefore permits one to understand the physical meaning of mass and heat transfer phenomena and to obtain a more accurate estimation of the related parameters.In an analogous way, the ENM has been used in order to fit the observed BTCs obtained from previous experiments on mass transport.
Section 3 shows the thermal tracer tests carried out on an artificially created fractured rock sample that has been used in previous studies to analyze nonlinear flow and non-Fickian transport dynamics in fractured formations (Cherubini et al., 2012(Cherubini et al., , 2013a(Cherubini et al., , b, c, 2014)).
In Sect. 4 have been reported the interpretation of flow and transport experiments together with the fitting of BTCs and interpretation of estimated model parameters.In particular, the obtained thermal BTCs show more enhanced early arrival and long tailing than solute BTCs.
The travel time for solute transport is an order of magnitude lower than for heat transport experiments.Thermal convective velocity is thus more delayed with respect to solute transport.The thermal dispersion mechanism dominates heat propagation in the fractured medium in the carried out experiments and thus cannot be neglected.
For mass transport the presence of the secondary path and the nonlinear flow regime are the main factors affecting non-Fickian behavior observed in experimental BTCs, whereas for heat transport the non-Fickian nature of the experimental BTCs is governed mainly by the heat exchange mechanism between the fracture network and the surrounding matrix.The presence of a nonlinear flow regime gives rise to a weak growth on heat transfer phenomena.
Section 5 reports some practical applications of the knowledge acquired from this study on the convective heat transport in fractured media for exploiting heat recovery and heat dissipation.Furthermore the estimation of the average effective thermal conductivity suggests that there is a solid thermal resistance in the fluid-to-solid heat transfer processes due to the rock-fracture size ratio.This result matches previous analyses (Pastore et al., 2015) in which a lower heat dissipation with respect to Tang's solution in correspondence to the single fracture surrounded by a matrix with more limited heat capacity has been found.With few exceptions, any fracture can be envisioned as two rough surfaces in contact.In cross section the solid areas representing asperities might be considered the grains of porous media.
Therefore, in most studies examining hydrodynamic processes in fractured media, the general equations describing flow and transport in porous media are applied, such as Darcy's law, which depicts a linear relationship between the pressure gradient and fluid velocity (Whitaker, 1986;Cherubini and Pastore, 2010).
However, this linearity has been demonstrated to be valid in low flow regimes (Re < 1).For Re > 1 a nonlinear flow behavior is likely to occur (Cherubini et al., 2013d).
When Re 1, a strong inertial regime develops that can be described by the Forchheimer equation (Forchheimer, 1901): where x (m) is the coordinate parallel to the axis of the single fracture (SF), p (ML ) is the convective velocity, ρ (ML −3 ) is the density and β (L −1 ) is called the inertial resistance coefficient, or non-Darcy coefficient.
It is possible to express the Forchheimer law in terms of hydraulic head h (L): The coefficients a (TL −1 ) and b (TL −2 ) represent the linear and inertial coefficient, respectively, equal to The relationship between hydraulic head gradient and flow rate Q (L 3 T −1 ) can be written as The coefficients a (TL −3 ) and b (T 2 L −6 ) can be related to a and b : where ω eq (L 2 ) is the equivalent cross-sectional area of SF.

Heat transfer by water flow in single fractures
Fluid flow and heat transfer in a single fracture (SF) undergo advective, diffusive and dispersive phenomena.Dispersion is caused by small-scale fracture aperture variations.Flow channeling is one example of macrodispersion caused by preferred flow paths, in that mass and heat tend to migrate through the portions of a fracture with the largest apertures.
In fractured media another process is represented by diffusion into the surrounding rock matrix.Matrix diffusion attenuates the mass and heat propagation in the fractures.
According to the boundary layer theory (Fahien, 1983), solute mass transfer q M (ML −2 ) per unit area at the fracturematrix interface (Wu et al., 2010) is given by where c f (ML −3 ) is the concentration across fractures, c m (ML −3 ) is the concentration of the matrix block surfaces, D m (LT −2 ) is the molecular diffusion coefficient, and δ (m) is the thickness of the boundary layer (Wu et al., 2010).For small fractures, δ may become the aperture w f (m) of the SF.
In an analogous manner, the specific heat transfer flux q H (MT −3 ) at the fracture-matrix interface is given by where T f (K) is the temperature across fractures, T m (K) is the temperature of the matrix block surfaces, and k m (MLT −3 K −1 ) is the thermal conductivity.
The continuity conditions at the fracture-matrix interface require a balance between mass transfer rate and mass diffused into the matrix described as where z (m) is the coordinate perpendicular to the fracture axis and w f is the aperture of the fracture.
In the same way, the specific heat flux must be balanced by heat diffused into the matrix described as The effective diffusion coefficient takes into account the fact that diffusion can only take place through pore and fracture openings because mineral grains block many of the possible pathways.The effective thermal conductivity of a formation consisting of multiple components depends on the geometrical configuration of the components as well as on the thermal conductivity of each.
The effective terms (D e instead of D m and k e instead of k m ) have been introduced in order to include the effect of various system parameters such as fluid velocity, porosity, surface area, and roughness that may enhance the mass and heat transfer effect.For instance, when large flow velocity occurs, convective transport is stronger along the center of the fracture, enhancing the concentration or temperature gradient at the fracture-matrix interface.As is known, roughness plays an important role in increasing mass or heat transfer because of increasing turbulent flow conditions.
According to Bodin (2007) the governing equation for the 1-D advective-dispersive transport along the axis of a semiinfinite fracture with 1-D diffusion in the rock matrix, in perpendicular direction to the axis of the fracture, is where D f (L 2 T −1 ) is the dispersion in the fracture.The latter mainly depends on two processes: Aris-Taylor dispersion and geometrical dispersion.Previous experiments (Cherubini et al., 2012a(Cherubini et al., , b, c, 2014) ) show that, due to the complex geometrical and topological characteristics of the fracture network that create tortuous flow paths, Aris-Taylor dispersion may not develop.A linear relationship has been found between velocity and dispersion, so geometrical dispersion is mostly responsible for the mixing process along the fracture: where α LM (L) is the dispersivity coefficient for mass transport.
Assuming that fluid flow velocity in the surrounding rock matrix is equal to zero, the equation for the conservation of heat in the matrix is given by where D a is the apparent diffusion coefficient of the solute in the matrix expressed as a function of the matrix porosity θ m , D a = D e /θ m (Bodin et al., 2007).Tang et al. (1981) presented an analytical solution for solute transport in a semi-infinite single fracture embedded in a porous rock matrix with a constant concentration at the fracture inlet (x = 0) equal to c 0 (ML −3 ) and with an initial concentration equal to zero.The solute concentration in the fracture c f and in the matrix c m are as follows: where s is the integral variable of the Laplace transform and L (L) is the length of SF; the v, A, β 2 and B coefficients are expressed as follows: whereas the gradient of c m at the interface z Defining the residence time as the average amount of time that the solute spends in the system, on the basis of these analytical solutions the probability density function (PDF) of the solute residence time in the single fracture in the Laplace space can be expressed as Assuming that density and heat capacity are constant in time, the heat transport conservation equation in SF can be expressed as follows: where ρ w (ML −3 ) and C w (L 2 T 2 K −1 ) represent the density and the specific heat capacity of the fluid within SF, respectively.D f for heat transport assumes the following expression: where λ L is the thermodynamic dispersion coefficient (MLT −3 K −1 ).Sauty et al. (1982) and de Marsily (1986) proposed an expression for the thermal dispersion coefficient where the thermal dispersion term varies linearly with velocity and depends on the heterogeneity of the medium, as for solute transport: where k 0 is the bulk thermal conductivity (MLT −3 K −1 ) and α LH (L) is the longitudinal thermal dispersivity.
The heat transport conservation equation in the matrix is expressed as follows: Note that the governing equations of heat and mass transport highlight similarities between the two processes; thus, Tang's solution can also be used for heat transport.
Three characteristic timescales can be defined: where L (L) is the characteristic length; t u (T), t d (T) and t e (T) represent the characteristics timescales of convective transport, dispersive transport and loss of the mass or heat into the surrounding matrix.
The relative effect of dispersion, convection and matrix diffusion on mass or heat propagation in the fracture can be evaluated by comparing the corresponding timescale.
Peclet number P e is defined as the ratio between dispersive (t d ) and convective (t u ) transport times: At high Peclet numbers transport processes are mainly governed by convection, whereas at low Peclet numbers it is mainly dispersion that dominates.Another useful dimensionless number, generally applied in chemical engineering, is the Damköhler number that can be used in order to evaluate the influence of matrix diffusion on convection phenomena.Da relates the convection timescale to the exchange timescale: where α (T −1 ) is the exchange rate coefficient corresponding to Note that the inverse of t e has the same meaning as the exchange rate coefficient α (T −1 ).
When t e values are of the same order of magnitude as the transport time t u (Da ∼ = 1), diffusive processes in the matrix are more relevant.In this case concentration or temperature distribution profiles are 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.By contrast, when t e t u (Da 1), the fracture-matrix exchange is rapid, there is instantaneous equilibrium between the fracture and the matrix, and they have the same concentration or temperature.These two circumstances close the standard advective-dispersive transport equation.
The product between Pe and Da represents another dimensionless group which is a measure of transport processes: When P e × Da increases, t e decreases more rapidly than t d , and subsequently the mass or heat diffusion into the matrix may be dominant on the longitudinal dispersion.

Explicit Network Model (ENM)
The 2-D Explicit Network Model (ENM) depicts the fractures as 1-D pipe elements forming a 2-D pipe network, and therefore expressly takes the fracture network geometry into account.The ENM permits one to understand the physical meaning of flow and transport phenomena and therefore to obtain a more accurate estimation of flow and transport parameters.
With the assumption that a j th SF can be schematized by a 1-D pipe element, the Forchheimer model can be used to write the relationship between head loss h j (L) and flow rate Q j (L 3 T −1 ) in finite terms: where L j (L) is the length of j th SF, and a (TL −3 ) and b (T 2 L −6 ) represent the Forchheimer parameters written in finite terms.The term in the square brackets constitutes the resistance to flow R j Q j (TL −2 ) of j th SF.
In case of steady-state conditions and for a simple 2-D fracture network geometry, a straightforward manner can be applied to obtain the solution of a flow field by applying the first and second Kirchhoff laws.
In a 2-D fracture network, fractures can be arranged in series and/or in parallel.Specifically, in a network in which fractures are set in a chain, the total resistance to flow is calculated by simply adding up the resistance values of each single fracture.The flow in a parallel fracture network breaks up, with some flowing along each parallel branch and recombining when the branches meet again.In order to estimate the total resistance to flow, the reciprocals of the resistance values have to be added up, and then the reciprocal of the total has to be calculated.The flow rate Q j across the generic fracture j of the parallel network can be calculated as (Cherubini et al., 2014) where ) is the sum of the mass flow rates at fracture intersections in correspondence to the inlet bond of j fracture, whereas the term in square brackets represents the probability of water distribution of j fracture P Q,j .
Once the flow field in the fracture network is known, to obtain the PDF at a generic node, the PDFs of each elementary path that reaches the node have to be summed up.They can be calculated as the convolution product of the PDFs of each single fracture composing the elementary path.
Definitely, the BTC describing the concentration in the fracture as a 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 and c inj (ML −3 ) is the concentration injection function, * is the convolution operator, L −1 represents the inverse Laplace transform operator, N p is the number of the paths reaching the node, n f,i is the number of the SF belonging to the elementary path ith, and P M,j and (s) are the mass distribution probability and the PDF in the Laplace space of the generic j th SF, respectively.The inverse Laplace transform L −1 can be solved numerically using the Abate and Ward (2006) algorithm.
In the same way the BTC T f which describes the temperature in the fracture as a function of time at the generic node can be written as where T 0 (K) is the initial temperature, T inj (K) is the temperature injection function and P H,j is the heat distribution probability.P M,j and P H,j can be estimated as the probabilities of the mass and heat distribution at the inlet bond of each individual SF, respectively.The mass and heat distribution is proportional to the correspondent flow rates: Note that if Eq. ( 38) is valid, the probability of water distribution is equal to the probabilities of mass and heat distribution (term in square brackets in Eq. 34).Therefore, the ENM model regarding each SF can be described by four parameters (u f,j , D f,j , α j , P Q,j ).
3 Material and methods

Description of the experimental apparatus
The heat transfer tests have been carried out on the experimental apparatus previously employed to perform flow and tracer transport experiments at bench scale (Cherubini et al. 2012(Cherubini et al. , 2013a(Cherubini et al. , b, c, 2014)).However, the apparatus has been modified in order to analyze heat transport dynamics.Two thermocouples have been placed at the inlet and the outlet of a selected fracture path of the limestone block with parallelepiped shape (0.6 × 0.4 × 0.08 m 3 ) described in previous studies.A TC-08 Thermocouple Data Logger (pico Technology) with a sampling rate of 1 s has been connected to the thermocouples.An extruded polystyrene panel with thermal conductivity equal to 0.034 Wm −1 K −1 and thickness 0.05 m has been used to thermally insulate the limestone block which has then been connected to a hydraulic circuit.The head loss between the upstream tank connected to the inlet port and the downstream tank connected to the outlet port drives flow of water through the fractured block.
An ultrasonic velocimeter (DOP3000 by Signal Processing) has been adopted to measure the instantaneous flow rate that flows across the block.An electric boiler with a volume of 10 −2 m 3 has been used to heat the water.In a flow cell located in correspondence to the outlet port, a multiparametric probe is positioned for the instantaneous measurement of pressure (dbar), temperature ( • C) and electric conductivity (µS cm −1 ). Figure 1a shows the fractured block sealed with epoxy resin and Figure 1b shows the thermal insulated fractured block connected to the hydraulic circuit, whereas the schematic diagram of the experimental apparatus is shown in Fig. 2.

Flow experiments
The average flow rate through the selected path can be evaluated as where S 1 (L 2 ) is the cross-sectional area of the flow cell, and t = t 1 − t 0 is the time for the flow cell to be filled from h 0 (L) and h 1 (L).To calculate the head loss between the upstream tank and the flow cell, the following expression is adopted: where h c is the hydraulic head measured in the upstream tank.Several tests have been carried out varying the control head, and in correspondence to each value of the average flow rate and head loss, the average resistance to flow has been determined as As an initial condition, a specific value of the hydraulic head difference between the upstream tank and the downstream tank has been assigned.At t = 0, valve a is closed so that the hydrostatic head inside the block assumes the same value as the one in the downstream tank.At t = 10 s, valve a is opened.

Solute and temperature tracer tests
For the solute tracer test at time t = 60 s by means of a syringe, a mass of 5 × 10 −4 kg sodium chloride is injected into the inlet port.Due to the very short source release time, the instantaneous source assumption can be adopted which assumes the source of a solute as an instantaneous injection (pulse).The multiparametric probe located within the flow cell measures the solute BTC.
As concerns thermal tracer tests at the time t = 60 s, the valve d is opened, while the valve c is closed.In such a way a step temperature function in correspondence to the inlet port T inj (t) is imposed and measured by the first thermocouple.The other thermocouple located inside the outlet port is used to measure the thermal BTC.
The ultrasonic velocimeter is used in order to measure the instantaneous flow rate, whereas a multiparametric probe located at the outlet port measures the pressure and the electric conductivity.

Flow characteristics
The Kirchhoff laws have been used in order to estimate the flow rates flowing in each single fracture.In Fig. 3 a sketch of the 2-D pipe conceptualization of the fracture network is reported.
The resistance to flow of each SF can be evaluated as the square bracket in Eq. ( 34).For simplicity the linear and nonlinear terms have been considered constant and equal for each SF.
The resistance to flow for the whole fracture network R Q can be evaluated as the sum of the resistance to flow of each SF arranged in a chain and the total resistance of the parallel branches equal to the reciprocal of the sum of the reciprocal of the resistance to flow of each parallel branch: where R j with j = 1-9 represents the resistance to flow of each SF, Q 0 is the injection flow rate, and Q 1 and Q 2 are the flow rates flowing in parallel branches 6 and 3-5, respectively.
The flow rate Q 1 is determined in an iterative manner using the following iterative equation derived by Eq. ( 35) at node 3: whereas the flow rate Q 2 is determined merely as The linear and nonlinear terms representative of the whole fracture network have been estimated by matching the average experimental resistance to flow resulting from Eq. ( 41) with resistance to flow estimated from Eq. ( 42).
The linear and nonlinear terms are equal, respectively, to a = 7.345 × 10 4 sm −3 and b = 11.65 × 10 9 s 2 m −6 .Inertial forces dominate viscous ones when the Forchheimer number (Fo) is higher than one.Fo can be evaluated as the ratio between the nonlinear loss bQ 2 and the linear loss (aQ).The critical flow rate Q crit which represents the value of the  flow rate for which Fo = 1 is derived as the ratio between a and b resulting in Q crit = 6.30× 10 −6 m 3 s −1 .
Because of the nonlinearity of flow, varying the inlet flow rate Q 0 , the ratio between the flow rates Q 1 and Q 2 flowing, respectively, in branches 6 and 3-5 is not constant.When Q 0 increases, Q 2 increases faster than Q 1 .The probability of the water distribution of the branch 6 P Q,6 is evaluated as the ratio between Q 0 and Q 1 , whereas the probability of the water distribution of branches 3-5 is equal to P Q,3−5 = 1 − P Q,6 .

Fitting of breakthrough curves and interpretation of estimated model parameters
The behavior of mass and heat transport has been compared by varying the injection flow rates.In particular, 21 tests in the range 1.83 × 10 −6 -1.26 × 10 −5 m 3 s −1 (Re in the range 17.5-78.71)for heat transport have been performed and compared with the 55 tests in the range 1.32-8.34× 10 −6 m 3 s −1 (Re in the range 8.2-52.1)for solute transport presented in previous studies.
The observed heat and mass BTCs for different flow rates have been individually fitted using the ENM approach presented in Sect.2.3.For simplicity, the transport parameters u f , D f and α are assumed equal for all branches of the fracture network.The probability of mass and heat distribution is assumed equal to the probability of water distribution.
The experimental BTCs are fitted using Eqs.( 36) and (37) for mass and heat transport, respectively.Note that, for mass transport, c inj (t) supposing the instantaneous injection condition becomes a Dirac delta function.
The determination coefficient (r 2 ) and the root mean square error (RMSE) have been used in order to evaluate the goodness of fit.
Tables 1 and 2 show the values of transport parameters, the RMSE and r 2 for mass and heat transport, respectively.Furthermore Figs. 4 and 5 show the fitting results of BTCs for some values of Q 0 .
The results presented in Tables 1 and 2 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.Regarding the transfer rate coefficient α, it assumes very low values for mass transport relative to the convective velocity.Instead, for heat transport the exchange rate coefficient is on the same order of magnitude of the convective velocity and, considering a characteristic length equal to L = 0.601 m corresponding to the length of the main path of the fracture network, the effect of dual porosity is very strong and cannot be neglected relative to the investigated injection flow range.Both mass and heat transport show a satisfactory fitting.In a particular manner, RMSE varies in the range 0.0015-0.0180for mass transport and in the range 0.0030-0.236for heat transport, whereas r 2 varies in the range 0.9863-0.9987for mass transport and in the range 0.0963-0.9998for heat transport.
In order to investigate the different behavior between mass and heat transport, the relationships between injection flow rate and the transport parameters have been analyzed.In Fig. 6 the relationship between u f and Q 0 is reported, whereas in Figs.7 and 8 the dispersion coefficient D f and the exchange term α as a function of u f , respectively, are reported.The figures show a very different behavior between mass and heat transport.
Regarding mass transport experiments according to previous studies (Cherubini at al., 2013a(Cherubini at al., , b, c, 2014)), Fig. 5 shows that for values of Q 0 higher than 4 × 10 −6 m 3 s −1 u f increases less rapidly.This behavior was due to the presence of inertial forces that gave rise to a retardation effect on solute transport.
Instead, Fig. 7 shows a linear relationship between u f and D f , suggesting that inertial forces did not exert any effect on dispersion and that geometrical dispersion dominates the Aris-Taylor dispersion.
In the same way as for mass transport, for heat transfer a linear relationship is evident between dispersion and convective velocity.Even if heat convective velocity is lower than solute advective velocity, the longitudinal thermal dispersivity assumes higher values than the longitudinal solute dispersivity.Also, for heat transport experiments, a linear relationship between u f and D f has been found.
Figure 8 shows the exchange rate coefficient α as a function of the convective velocity u f for both mass and heat transport.
Regarding the mass transport, the estimated exchange rate coefficient α is much lower than the convective velocity.These results suggest that in the case study fracture-matrix Figure 5. Fitting of BTCs at different injection flow rates using the ENM with Tang's solution for heat transport.The blue curve is the temperature observed at the inlet port used as the temperature injection function, the red square curve is the observed temperature at the outlet port, and the black continuous curve is the simulated temperature at the outlet port.
exchange is very slow and that it may not influence mass transport.Non-Fickian behavior observed in the experimental BTCs is therefore dominated mainly by the presence of inertial forces and the parallel branches.
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 fracture-matrix exchange phenomena resulting in a significant retardation effect.Once the model parameters for each flow rate have been determined, the unit response function (f URF ), corresponding to the PDF obtained from impulsive injection of both solute and temperature tracers, is obtained.The unit response function can be characterized using the time moments and tail character analysis.
The mean residence time t m assumes the following expression: whereas the nth normalized central moment of distribution of f URF versus time can be written as The second moment µ 2 can be used in order to evaluate the dispersion relative to t m , whereas the skewness is a measure www.nonlin-processes-geophys.net/24/23/2017/ Nonlin.Processes Geophys., 24, 23-42, 2017     of the degree of asymmetry and is defined as follows: The tailing character t c can be described as where t fall denotes the duration of the falling limb defined as the time interval from the peak to the tail cutoff, which is the time when the falling limb first reaches a value that is 0.05 times the peak value.t rise is defined as the time interval from the first arrival to the peak.This quantity provides a measure of the asymmetry between the rising and falling  limbs.A value of t c significantly higher than 1 indicates an elongated tail compared to the rising limb (Cherubini et al., 2010).In Fig. 9 is reported the residence time versus the injection flow rates.The figure highlights that t m for heat transport is about 3 times higher than for mass transport.In a particular way, t m varies in the range 40.3-237.1 s for mass transport and in the range 147.8-506.9s for heat transport.This result still highlights that heat transport is more delayed than mass transport.
In the same way the skewness S (Fig. 10) and the tailing character t c (Fig. 11) are reported as a function of Q 0 .
A different behavior for heat and mass transport is observed for the skewness coefficient.For heat transfer the skewness shows a growth trend which seems to decrease af-   ter Q 0 = 3 × 10 −6 m 3 s −1 .Its mean value is equal to 2.714.
For solute transport the S does not show a trend, and assumes a mean value equal to 2.018.The tailing character does not exhibit a trend for either mass and heat transport.In either cases t c is significantly higher than 1, specifically 7.70 and 30.99 for mass and heat transport, respectively.
In order to explain the transport dynamics, the trends of dimensionless numbers Pe and Da varying the injection flow rate have been investigated.Figure 12 shows the Pe as a function of Q 0 for both mass and heat experiments.As concerns mass experiments, Pe increases as Q 0 increases, assuming a constant value for high values (Pe = 7.5) of Q 0 .For heat transport a different behavior is observed, P e showing a con-  stant trend and being always lower than one.Even if the injection flow rate is relatively high, thermal dispersion is the dominating mechanism in heat transfer.
Figure 12 reports Da as a function of Q 0 .For mass transport Da assumes very low values, on the order of magnitude of 10 −4 .
The convective transport scale is very low with respect to the exchange transport scale; thus, the mass transport in each single fracture can be represented with the classical advection dispersion model.
As regards heat transport, Da assumes values around the unit showing a downward trend as injection flow rate increases, switching from higher to lower values than the unit.As injection flow rate increases, the convective transport timescale reduces more rapidly than the exchange timescale.
These arguments can be explained because the relationships between Q 0 and u f show a change in slope when Da becomes lower than the unit.In other words, when Da is higher than the unit, the exchange between fracture and matrix dominates on the convective transport, giving rise to a more enhanced delay on heat transport; conversely, when Da is lower than one, convective transport dominates on fracture-matrix interactions and the delay effect is reduced.
Furthermore this effect is evident also in the trend observed in the graph S − Q 0 (Fig. 10).For values of Da lower than the unit, a change in slope is evident; the skewness coefficient increases more slowly.Thus for Da > 1 the early arrival and the tail effect of BTC increase more rapidly than for Da < 1.
Note that 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, which however assumes values lower than the unit.
Figure 14 shows the dimensionless group Pe × Da varying the injection flow rate.Regarding mass transport, Pe × Da is on the order of magnitude of 10 −3 , confirming the fact that the fracture-matrix interaction can be neglected relative to the investigated range of injection flow rates.For heat transport, Pe × Da assumes values just below the unit, with a downward trend as Q 0 increases.t d and t e have the same order of magnitude.
In order to find the optimal conditions for heat transfer in the analyzed fractured medium, the thermal power exchanged per unit temperature difference Q (ML 2 T −1 K −1 ) for each injection flow rate in quasi-steady-state conditions can be estimated.The thermal power exchanged can be writ-ten as The outlet temperature T out can be evaluated as a function of the f URF using the following expression: Substituting Eq. ( 50) into Eq.( 49), the thermal power exchanged per unit temperature difference is Figure 15 shows the similarities between the relationship Q/ T inj − T 0 − Q 0 (Fig. 15a) and Da − Q 0 (Fig. 14b).Higher Da values correspond to higher values of Q/ T inj − T 0 .The thermal power exchanged increases as the Damköhler number increases, as shown in Fig. 15c.These results highlight that for the observed case study the optimal condition for thermal exchange in the fractured medium is obtained when the exchange timescale is lower than the convective transport scale, or rather when the dynamics of fracture-matrix exchange are dominant on the convective ones.Moreover, in a similar way to Da, Q/ T inj − T 0 shows a weak growth trend when Q 0 exceeds Q crit .This means that the nonlinear flow regime improves the fracture-matrix thermal exchange; however, at high values of injection flow rates, convective and dispersion timescales are less than the exchange timescale.Nevertheless, these results have been observed in a small range of Da numbers close to the unit.In order to generalize these results, a larger range of Da numbers should be investigated.
In order to estimate the effective thermal conductivity coefficient k e , the principle of conservation of heat energy can be applied to the whole fractured medium.Neglecting the heat stored in the fractures, the difference between the heat measured at the inlet and at the outlet must be equal to the heat diffused into the matrix: where A f is the whole surface area of the whole active fracture network and the gradient of T m can be evaluated according to Eq. ( 19) using temperature instead of concentration as a variable.Then the average effective thermal conductivity k e can be obtained as  The average effective thermal conductivity has been estimated for each injection flow rate (Fig. 16) and assumes a mean value equal to k e = 0.1183 Wm −1 K −1 .The estimated k e is 1 order of magnitude lower than the thermal conductivity coefficient reported in the literature (Robertson, 1988).
Fractured media have a lower capacity for diffusion, as opposed to Tang's model, which has unlimited capacity.There is a solid thermal resistance in the fluid to solid heat transfer processes which depends on the rock-fracture size ratio.This result is coherent with previous analyses on heat transfer carried out on the same rock sample (Pastore et al., 2015).In this study Pastore et al. (2015) found that the ENM model failed to model the behavior of heat transport, in correspondence to parallel branches where the hypothesis of Tang's solution of a single fracture embedded in a porous medium having unlimited capacity cannot be considered valid.In parallel branches the observed BTCs are char-acterized by less retardation of heat propagation as opposed to the simulated BTCs.

Conclusions
Aquifers offer a possibility of exploiting geothermal energy by withdrawing the heat from groundwater by means of a heat pump and subsequently supplying the water back into the aquifer through an injection well.In order to optimize the efficiency of the heat transfer system and minimize the environmental impacts, it is necessary to study the behavior of convective heat transport especially in fractured media, where flow and heat transport processes are not well known.
Laboratory experiments on the observation of mass and heat transport in a fractured rock sample have been carried out in order to analyze the contribution of thermal dispersion in heat propagation processes, the contribution of nonlinear flow dynamics to the enhancement of thermal matrix diffu- sion, and finally the optimal heat recovery and heat dissipation strategies.
The parameters that control mass and heat transport have been estimated using the ENM model based on Tang's solution.
Heat transport shows a very different behavior compared to mass transport.The estimated transport parameters show differences of several orders of magnitude.Convective thermal velocity is lower than solute velocity, whereas thermal dispersion is higher than solute dispersion, mass transfer rate assumes a very low value, suggesting that fracture-matrix mass exchange can be neglected.Non-Fickian behavior of observed solute BTCs is mainly due to the presence of the secondary path and the nonlinear flow regime.Contrarily, heat transfer rate is comparable with convective thermal velocity giving rise to a retardation effect on heat propagation in the fracture network.
The discrepancies detected in transport parameters are moreover observable through the time moment and tail character analysis which demonstrate that the dual porosity behavior is more evident in the thermal BTCs than in the solute BTCs.
The dimensionless analysis carried out on the transport parameters proves that, as the injection flow rate increases, thermal convection timescale decreases more rapidly than the thermal exchange timescale, explaining the reason why the relationship Q 0 − u f shows a change in slope for Da lower than the unit.
Thermal dispersion dominates heat transport dynamics and the Peclet number, and the product between the Peclet number and the Damköhler number is almost always less than the unit.
The optimal conditions for thermal exchange in a fracture network have been investigated.The power exchanged increases in a potential way as Da increases in the observed range.
The Explicit Network Model is an efficient computation methodology to represent flow, mass and heat transport in fractured media, as 2-D and/or 3-D problems are reduced to resolving a network of 1-D pipe elements.Unfortunately, in field case studies, it is difficult to obtain full knowledge of the geometry and parameters such as the orientations and aperture distributions of the fractures needed by the ENM, even by means of field investigation methods.However, in real case studies the ENM can be coupled with continuum models in order to represent greater discontinuities with respect to the scale of study, which generally gives rise to preferential pathways for flow, mass and heat transport.
A method to represent the topology of the fracture network is represented by multifractal analysis as discussed in Tijera at al. (2009) and Tarquis at al. (2014).
This study has permitted one to detect the key parameters to design devices for heat recovery and heat dissipation that exploit the convective heat transport in fractured media.
Heat storage and transfer in fractured geological systems is affected by the spatial layout of the discontinuities.
Specifically, the rock-fracture size ratio which determines the matrix block size is a crucial element in determining matrix diffusion on the fracture-matrix surface.
The estimation of the average effective thermal conductivity coefficient shows that it is not efficient to store thermal energy in rocks with high fracture density because the fractures are surrounded by a matrix with more limited capacity for diffusion giving rise to an increase in solid thermal resistance.In fact, if the fractures in the reservoir have a high density and are well connected, such that the matrix blocks are small, the optimal conditions for thermal exchange are not reached, as the matrix blocks have a limited capability to store heat.
On the other hand, isolated permeable fractures will tend to lead to more distribution of heat throughout the matrix.
Therefore, subsurface reservoir formations with large porous matrix blocks will be the optimal geological formations to be exploited for geothermal power development.
The study could help to improve the efficiency and optimization of industrial and environmental systems, and may provide a better understanding of geological processes involving transient heat transfer in the subsurface.
Future developments of the current study will be carrying out investigations and experiments aimed at further deepening of the quantitative understanding of how fracture arrangement and matrix interactions affect the efficiency of storing and dissipating thermal energy in aquifers.This could be achieved by means of using different formations with different fracture density and matrix porosity.

C
.Cherubini et al.:  Laboratory experimental investigation of heat transport in fractured Figure 1.(a) Fractured block sealed with epoxy resin.(b) Thermal insulated fracture block connected to the hydraulic circuit.

Figure 2 .
Figure 2. Schematic diagram of the experimental setup.

Figure 3 .
Figure 3. Two-dimensional pipe network conceptualization of the fracture network of the fractured rock block in Fig. 1.Q 0 is the injection flow rate; Q 1 and Q 2 are the flow rates that flow in parallel branches 6 and 3-5, respectively.

Figure 4 .
Figure 4. Fitting of BTCs at different injection flow rates using the ENM with Tang's solution for mass transport.The green square curve is the observed specific mass flux at the outlet port; the continuous black line is the simulated specific mass flux.

Figure 6 .
Figure 6.Velocity u f (m s −1 ) as a function of the injection flow rate Q 0 (m 3 s −1 ) for the ENM with Tang's solution for both mass and heat transport.

Figure 7 .
Figure 7. Dispersion D f (m s −2 ) as a function of velocity u f (m s −1 ) for the ENM with Tang's solution for both mass and heat transport.

Figure 8 .
Figure 8. Transfer coefficient α (s −1 ) as a function of velocity u f (m s −1 ) for both mass and heat transport.

Figure 9 .
Figure 9. Mean travel time t m (s) as a function of the injection flow rate for both mass and heat transport.

Figure 10 .
Figure 10.Skewness as a function of the injection flow rate for both mass and heat transport.

Figure 11 .
Figure 11.Tailing character t c as a function of the injection flow rate for both mass and heat transport.

Figure 12 .
Figure 12.Peclet number as a function of the injection flow rate Q 0 (m 3 s −1 ) for both mass and heat transport.

Figure 13 .
Figure 13.Da number as a function of the injection flow rate Q 0 (m 3 s −1 ) for both mass and heat transport.

Figure 14 .
Figure 14.Pe × Da number as a function of the injection flow rate Q 0 (m 3 s −1 ) for both mass and heat transport.

Figure 15 .
Figure 15.Heat power exchanged per difference temperature unit Q/(T inj − T 0 ) as a function of the injection flow rate Q 0 (m 3 s −1 ) (a), Damköhler number Da as a function of the injection flow rate (b), and power exchanged per difference temperature unit as a function of the Damköhler number (c).

Table 1 .
Estimated values of parameters, RMSE, and determination coefficient r 2 for the ENM with Tang's solution at different injection flow rates for mass transport.

Table 2 .
Estimated values of parameters, RMSE, and determination coefficient r 2 for the ENM with Tang's solution at different injection flow rates for heat transport.