Thermal damping and retardation in karst conduits

Water temperature is a non-conservative tracer in the environment. Variations in recharge temperature are damped and retarded as water moves through an aquifer due to heat exchange between water and rock. However, within karst aquifers, seasonal and short-term fluctuations in recharge temperature are often transmitted over long distances before they are fully damped. Using analytical solutions and numerical simulations, we develop relationships that describe the effect of flow path properties, flow-through time, recharge characteristics, and water and rock physical properties on the damping and retardation of thermal peaks/troughs in karst conduits. Using these relationships, one can estimate the thermal retardation and damping that would occur under given conditions with a given conduit geometry. Ultimately, these relationships can be used with thermal damping and retardation field data to estimate parameters such as conduit diameter. We also examine sets of numerical simulations where we relax some of the assumptions used to develop these relationships, testing the effects of variable diameter, variable velocity, open channels, and recharge shape on thermal damping and retardation to provide some constraints on uncertainty. Finally, we discuss a multitracer experiment that provides some field confirmation of our relationships. High temporal resolution water temperature data are required to obtain sufficient constraints on the magnitude and timing of thermal peaks and troughs in order to take full advantage of water temperature as a tracer.


Introduction
Much of the flow and transport through karst aquifers occurs via conduits (Atkinson, 1977b;Worthington, 1999;Worthington et al., 2000).These preferential flow paths occur in all karst aquifers, but most are poorly characterized or unknown.Hydrogeological investigations of karst aquifers frequently employ hydrographs, chemographs, and thermographs collected from boreholes and springs.Ideally, these data could be used to provide flow path information, and this characterization would facilitate models that are more capable of predicting flow and transport through these systems.
Hydrographs have been analyzed for more than a century (Boussinesq, 1903(Boussinesq, , 1904;;Maillet, 1905) to characterize flow recession, determine aquifer characteristics, or predict discharge with time (e.g., see summaries in Hall, 1968;Tallaksen, 1995;Jeannin and Sauter, 1998;Dewandel et al., 2003;Ford and Williams, 2007).However, hydrographs provide minimal information about conduit geometry (Covington et al., 2009), and interpretations of karst aquifer structures based on hydrograph analysis are problematic because of the relatively strong control of rainfall frequency on hydrograph shape (Jeannin and Sauter, 1998).Variations in specific conductance often occur with changes in localized recharge (Jakucs, 1959;Newson, 1971;Ternan, 1972;Atkinson, 1977a, b;Worthington et al., 1992;White, 2002).While chemical modification of electrical conductivity signals due to dissolving calcite could theoretically be used to constrain the geometry of flow paths with hydraulic diameters on the mm to cm scale, electrical conductivity provides little in-formation about conduits with diameters on the meter scale and larger because these larger flow paths produce negligible chemical modification of localized recharge from dissolution (Covington et al., 2012).
Conduits facilitate fast flow-through times and may enable thermal perturbations to reach a spring (e.g., Benderitter et al., 1993;Bundschuh, 1997;Martin and Dean, 1999;Screaton et al., 2004;Luhmann et al., 2011).These perturbations are modified as water flows through the system, and the modification is sensitive to conduit geometry (Renner, 1996;Liedl et al., 1998;Liedl and Sauter, 1998).The modification occurs because of the heat exchange between water and rock, causing both damping (i.e., decrease in signal amplitude) and retardation (i.e., time lag of the signal) of recharge (Luhmann et al., 2012).Studies have also demonstrated thermal damping and retardation in porous media (e.g., Molson et al., 1992;Palmer et al., 1992;Markle and Schincariol, 2007) and fractures (e.g., Molson et al., 2007).The nonconservative nature of water temperature, even within fairly large conduits, facilitates estimates of conduit size via an analysis of input and output thermographs (Covington et al., 2011(Covington et al., , 2012;;Luhmann et al., 2012;Birk et al., 2014).Unlike chemical modification, the degree of thermal modification depends on the timescale of recharge variations.Shorter, storm-event thermal perturbations provide maximum information about conduits with hydraulic diameters on the meter scale; longer, seasonal thermal perturbations probe smaller, mm to cm scale flow paths (Covington et al., 2012).Previous work has also demonstrated that groundwater input into surface streams in karst terrains modifies the relationships between air and water temperatures (O' Driscoll and DeWalle, 2006).The extent of this modification will depend upon whether groundwater has had sufficient residence time to reach thermal equilibration (Luhmann et al., 2011;Covington et al., 2012).
In addition to correlations between thermal signals and conduit geometry, temperature peaks have been used as a simple and inexpensive means of estimating residence times within karst conduit systems when the timing of changes in recharge temperature is known (Martin and Dean, 1999;Birk et al., 2004;Screaton et al., 2004;Covington et al., 2011;Gunn, 2015).However, since heat exchange within a karst conduit introduces a retardation in the timing of the peak, residence times estimated using temperature will typically be longer than the true residence time.The magnitude of this error, and its functional relationship with conduit geometry and boundary conditions, have not been previously quantified, though Birk et al. (2004) noted that estimates of conduit volume based on temperature lags displayed significantly more scatter than estimates using electrical conductivity lags, and concluded that electrical conductivity provided a more reliable means of estimating travel times.
Recent work used temperature to identify water sources by employing a two-component mixing model (Doucette and Peterson, 2014).However, since heat exchange within a karst conduit dampens all thermal perturbations, there will be error in estimates of different water source fractions derived from models that assume conservative end-member temperatures.Temperature mixing models will typically overestimate contributions from background temperature sources and underestimate source waters that provide the thermal perturbations at the thermal peak/trough.Alternatively, during the thermal recession, the heated or cooled rock surrounding the conduit may potentially facilitate water temperatures that are no longer within the temperature range of the different water sources.The magnitude of these errors will depend upon the extent of water temperature change that occurs along the flow paths.
Our primary objective in this study is to demonstrate the effect of conduit geometry on thermal damping and retardation in karst conduits using both analytical solutions and numerical simulations.We also consider the effects of fluid flow velocity, recharge characteristics, and rock and water physical properties.A relationship between conduit geometry and thermal damping or retardation may ultimately be used to estimate conduit diameter given recharge temperature and down-gradient monitoring data that include water temperature and a conservative tracer.These relationships can also be used to estimate, and potentially correct for, errors in residence times or water source fractions derived from temperature pulses, and to understand how these errors vary with conduit properties and recharge.

Conceptual model
A simplified conceptual model of heat transport in a conduit through a karst aquifer is employed to provide a general understanding of thermal damping and retardation in karst conduits.Heat transport occurs via advection and dispersion in the conduit, conduction in the rock surrounding the conduit, and exchange between the conduit and rock.Both a circular conduit and a planar fracture or bedding plane are investigated.Velocity is imposed along the entire length of the conduit, and the following analysis generally assumes a number of constants to simplify the system, including water flow at a constant velocity in a conduit with a constant hydraulic diameter.However, the effects of velocity variations induced by recharge events and a conduit diameter that varies along its length are briefly considered.The analysis also generally employs conduits with full pipe flow, but open channel flow, where conduits are only partially full of water, is also considered in a few example cases.The conceptual model does not account for exchange of water between the conduit and matrix (or fractures or other conduits) or spatial velocity gradients or mobile/immobile regions within the conduit.Finally, the analysis employs a number of simple functions to approximate the shape of thermal perturbations produced in nature.Some limitations of our conceptual model are discussed in Sect.9.2.

Mathematical model
Temperature along a karst conduit as a function of time can be approximated by the 1-D heat advection-dispersion equation: where T w is the water temperature, t is time, D L is the longitudinal dispersivity, x is the longitudinal distance down the conduit, V is the water velocity, h conv is the convective heat transfer coefficient, ρ w is the density of water, c p,w is the specific heat of water at constant pressure, D H is the hydraulic diameter of the flow path, and T s is the conduit wall surface temperature.The terms on the right side of Eq. ( 1) describe heat dispersion, heat advection, and heat exchange with the surrounding rock.The convective heat transfer coefficient is given by where k w is the thermal conductivity of water and Nu is the dimensionless Nusselt number, which is the ratio of convective to purely conductive heat transfer through the convective boundary layer near the wall.Nu for turbulent flow is given by the empirically derived Gnielinski correlation (Incropera et al., 2007, Eq. 8.62): where f is the Darcy-Weisbach friction factor, Re = ρ w V D H /µ w is the dimensionless Reynolds number, Pr = c p,w µ w /k w is the dimensionless Prandtl number of water, and µ w is the dynamic viscosity of water.Conduction provides a strong control over heat exchange in karst conduits (Covington et al., 2011).Heat conduction in the rock surrounding a circular conduit with no energy generation can be described, using cylindrical symmetry, by the 2-D heat conduction equation where r is the radial distance from the conduit center, T r is the rock temperature, and α r = k r /(ρ r c p,r ) is the hereassumed isotropic and homogeneous thermal diffusivity of rock, with k r denoting the thermal conductivity, ρ r the density, and c p,r the specific heat.To represent heat transport in rock adjacent to a planar fracture or a bedding plane parting, we use translational symmetry, and the heat conduction equation becomes where y is the distance from the fracture center.Furthermore, heat exchange in a cylindrical conduit can be approximated by heat exchange in planar coordinates in many cases (Covington et al., 2011), permitting simpler planar simulations.The boundary conditions are and In planar coordinates, r in Eqs. ( 7) and ( 8) is replaced by y: and

Analytical solutions for damping and thermal retardation
Simplification of Eq. ( 1) allows derivation of several useful analytical solutions.Karst conduits are frequently advection dominated, with Peclet numbers of around 100 (Field and Nash, 1997).Therefore, neglecting longitudinal dispersivity will provide a reasonable approximation in many cases.This approximation will break down for particularly short duration pulses (Hauns et al., 2001), but is more likely to hold for the longer-term pulses typically found from natural perturbations.Neglecting longitudinal dispersivity results in For most relevant cases, where the timescale of the change in water temperature is not extremely short, the approximation h conv → ∞ is valid (Covington et al., 2011).In this case heat flow is limited by conduction in the rock, and one obtains a boundary condition T r (x, r or y = conduit wall, t) = T w (x, t). (12)

Sinusoidal solution for the planar case
Heat conduction in the rock along the length of the conduit (the x direction) is neglected, and thus, the equation for heat conduction in the rock becomes Equations ( 10) to ( 13) can be solved for the case of sinusoidally varying water temperature, allowing direct calculation of the thermal damping and retardation of the input wave.The retardation and damping produced for this sinusoidal upstream boundary condition provide significant insight into the response from many pulses found in natural settings, including, as will be seen later, a single isolated pulse.
For the sinusoidal solution, we employ a shifted temperature variable defined as where T r,∞ is the rock temperature at an infinite distance from the conduit axis.For an upstream boundary condition that is sinusoidal in time, the solution for the rock temperature is separable and has the functional form Since the water and rock temperatures are assumed equal at the boundary (Eq.12), Eq. ( 15) contains all of the temperature field.With a sinusoidal upstream boundary condition, the time varying component of the solution is also sinusoidal, T t (t) = T w,in e −iωt and T r (x, y, t) = T w,in X(x)Y (y)e −iωt , where T w,in is the peak amplitude of the input temperature variation.Using this in Eq. ( 13), combined with the boundary condition lim y→∞ T r = 0, leads to The function X(x) can then be derived using Eqs.( 10)-( 12), leading to where = ρ w c p,w /(ρ r c p,r ) is a ratio of the volumetric heat capacities of water and rock.This is an ordinary differential equation with constant coefficients and the solution is an exponential function X(x) = e −γ x , where For the water temperature at the conduit outlet, T w,out , this gives the solution where L is the conduit length.Since we are interested only in real solutions, we fix the phase and only look at the real part of the equation.
From this solution, one can directly derive both the retardation and damping experienced by each sinusoidal temperature peak.A peak in the output temperature is reached at a distance L (i.e., conduit length) downstream of the input at the time, t peak,out , when the imaginary part of the exponent is zero; that is, The fluid flow-through time through the conduit is t ft = L/V , and the retardation of the thermal peak, τ , is the difference As can been seen from the real part of γ , the damping of the upstream temperature peaks observed at the downstream end of the conduit (x = L) is given by This solution illustrates a thermal length scale, λ T,sin , that is appropriate for sinusoidal temperature variations in the input temperature, with λ T,sin is, to within a dimensionless factor of order 1, equivalent to the late time thermal length scale of Eq. ( 22) in Covington et al. (2012).

Sinusoidal solution for the cylindrical case
For heat conduction within the rock in the vicinity of a karst conduit with a cylindrical geometry, we again neglect conduction in the direction along the conduit (x) and, instead of Eq. ( 4), we use The solution remains separable, such that Again, we use sinusoidal T t (t) and get T r (x, r, t) = T w,in X(x)R(r)e −iωt .Substituting this into Eq.( 24) gives a Bessel equation whose solutions are Bessel functions.From the boundary condition lim r→∞ T r = 0 follows lim r→∞ R(r) = 0, which limits the solution space for R(r) to specific linear combinations of Bessel functions that are known as Hankel functions of the first kind, H i .The solution is As in the planar case, X(x) is obtained from Eq. ( 11) and has the form X(x) = e −γ x , where Because of the special functions, this solution is less useful analytically, but provides a straightforward means of calculating the output wave numerically.

Numerical integration of the planar case for arbitrary recharge temperature
As shown above, if the temperature at the input is T w,in (t) = e −iωt , then the temperature at the output is T w,out (t) = T w,in (t)X(L) = e −γ (ω)L−iωt .A general T w,in (t) can be expressed in terms of its Fourier transform, and the output temperature is then calculated as A Gaussian pulse is of particular interest, since this shape approximates many natural recharge events and is also the functional form we use for the simulations below.We use a Gaussian recharge function of the form where T w,in is T w at x = 0, T r,0 is the initial rock temperature (or T r,∞ ), R A is the recharge temperature amplitude, t peak,in is the peak time at x = 0, and σ controls the width of the thermal pulse.The Fourier transform is given by Therefore, the general solution for a Gaussian pulse can be calculated using the integral In practice, this equation, or Eq. ( 29) for the general case, must be integrated numerically.However, the Fourier transform solution provides an efficient means of numerically calculating thermographs.

Numerical simulations
In order to relax the somewhat restrictive assumptions required by the analytical solutions, and particularly to test the applicability of the sinusoidal analytical solutions to the propagation of isolated pulses, we present the results of numerical simulations of thermal pulses.These simulations solve the full version of Eqs. ( 1) and ( 4) or (5) for a variety of recharge and flow conditions, conduit geometries, thermal pulse shapes, rock and water physical properties, and also for open channel cases that include radiative heat exchange.For the majority of the simulation set, recharge temperature is varied according to the Gaussian function given in Eq. ( 30).
For each simulation, σ is defined to attain a desired recharge duration, R D , or full width at half maximum given by For the initial condition, T w and T r are set equal to T r,∞ or For most of the simulations, V is constant, although V varies between different simulations.f is approximated for most simulations using the von Kármán equation, where R = D H /2 is the conduit radius and is the roughness height (i.e., the average distance that irregularities on the rock wall protrude into the conduit).We set = 2.15 cm for all simulations.We also run simulations where f is calculated using the empirical Colebrook-White equation, and we find that simulation results are identical regardless of the equation used to determine f (Luhmann, 2011).
The COMSOL Multiphysics ® (Version 3.5) finite element package is used to solve the coupled heat advectiondispersion and conduction equations numerically.Using the Coefficient Form PDE mode in COMSOL, Eq. ( 1) is solved along a 1-D line, which represents a conduit (Fig. 1a) or fracture (Fig. 1b).Because of axial symmetry, a simulation of conduction in the rock surrounding a circular conduit with full pipe flow may be reduced to a 2-D axisymmetric problem.Thus, Eq. ( 4) is solved using COMSOL's Conduction Heat Transfer application mode with a 2-D axisymmetric rectangle for cylindrical simulations (Fig. 1a).Similarly, because of translational symmetry across the fracture plane, a simulation of conduction in the rock surrounding a water-filled fracture may be simplified to a 2-D planar problem.Thus, Eq. ( 5) is solved in a 2-D rectangle in Cartesian coordinates for planar simulations (Fig. 1b).The  7) or ( 9), respectively, and are set to background temperature.
Table 1.Default parameters used in simulations.

Parameter Value Unit
from the flow path wall, but the 2-D rectangle was generally discretized into 23 000 elements.COMSOL uses an implicit method to solve the system of equations.User-defined relative and absolute tolerances are compared to the estimated error to modify timestep duration to obtain the desired accuracy.The relative and absolute tolerances were set to ≤ 10 −6 and ≤ 10 −7 , respectively.Several example cases were run at higher spatial and temporal resolution, and produced the same results.
We conduct numerous simulations where we vary the parameters to consider their effect on thermal damping and retardation.Table 1 lists default values for parameters, but simulations were also run with other values, which are provided in the Supplement.The thermal transmission factor, F , which provides a means to quantify damping, is given by where T w,peak,out and T w,peak,in are the outlet and inlet peak water temperatures, respectively.The thermal peak retardation, τ , for each simulation is the time of peak temperature at the outlet minus the flow-through time (Eq.21).Though the notation here is in terms of temperature peaks, the same equations apply to temperature troughs.τ and F for all simulations are provided in the Supplement.

Thermal damping
The damping of temperature peaks in the simulations is dependent on the ratio L/λ T,sin .When the ratio L/λ T,sin is small, there is little damping of recharge signals.However, when the ratio L/λ T,sin is large, recharge signals undergo significant to complete damping.For the planar sinusoidal solution, the transmission factor, F , is given by Eq. ( 22).
In order to compare this analytical solution for damping of sinusoids with the simulations that contain Gaussian input thermographs, we need an approximate conversion between the angular frequency of the sinusoid, ω, and an appropriate analog for the Gaussian pulse.We use the relation ω = π/(C time R D ), which relates the period of the sinusoid to a multiple of the full width at half maximum of the Gaussian curve.The time conversion constant, C time , is treated as a fitting parameter.Using this approximation, and the definition of the transmission factor (Eq. 35), Eq. ( 22) can be rewritten

Big Cyl Small Cyl Slow Cyl Planar
Figure 2. A comparison of the transmission factors of peaks in the simulations of Gaussian temperature pulses against the modified form of the analytical solution for a sinusoidal input temperature (Eq.36).Cylindrical cases are corrected by an additional factor (Eq. 37) that is a function of the dimensionless parameter .These modified forms of the analytical solution provide a close fit to the simulation results for most cases.Big Cyl and Small Cyl indicate a conduit in cylindrical coordinates with a D H ≥ 1 m and a D H < 1 m, respectively.Slow Cyl indicates a conduit in cylindrical coordinates with a V ≤ 0.0352 m s −1 .Planar indicates a conduit in planar coordinates.The graph includes a 1 : 1 line. as For the planar simulations, and cylindrical simulations that are well approximated with the planar solution, we find that a value of C time ≈ 4 provides a tight fit to the transmission factors measured from the pulses in the simulations (Fig. 2).Covington et al. (2011) showed that the agreement between planar and cylindrical heat transport solutions was dependent on a dimensionless quantity, = (R 2 V )/(Lα r ), where V is a time-averaged or reference flow velocity, with cylindrical cases well approximated by the planar solution for 10.Similarly, here we find that Eq. ( 36) breaks down for cylindrical simulations with small .However, we also find that the error is strongly correlated with , and the damping in the cylindrical cases is well fit by a correction factor of the form This correction factor, with a value of C cyl ≈ 0.4, roughly accounts for the additional heat exchange.Though an approximation of the cylindrical solution given in Sect.4.2 might provide a more well-motivated correction, this equation produces acceptable results and is simpler to implement since it requires no calculation of Hankel functions.simulations and the values of F planar or F cyl for the planar and cylindrical simulations, respectively.

Thermal retardation
As for thermal damping, the thermal retardation of a Gaussian pulse can be approximated using the form of the sinusoidal solution along with a multiplicative time constant.For thermal retardation, we find that Eq. ( 21) provides a good approximation to the simulated cases (Fig. 3) with a choice of ω ≈ π/R D , such that While this relation provides an excellent fit to the planar cases, and most of the cylindrical cases, cylindrical cases with small values of do produce some scatter.This scatter is sufficiently small that we do not attempt to develop a correction for it.There is also some scatter associated with simulations with relatively slow velocities.This scatter is likely caused by numerical dispersion.

Relaxation of additional assumptions
Our analysis thus far, including the simulations, employs a number of simplifications or approximations, such as constant conduit diameter and constant flow velocity.Here, we run simulations that relax these assumptions to examine potential uncertainty in the F and τ relationships.We consider the effect of variable diameter or flow velocity within an individual conduit and also run open channel simulations, where a conduit is only partially filled with water and radiative heat exchange occurs.Finally, we consider other functions that approximate the shape of recharge thermographs in nature to examine whether different shapes produce significantly different values of damping or retardation.

Variable hydraulic diameter
The conduit hydraulic diameter, D H , typically changes along a karst flow path.If this occurs, the thermal signal at the monitoring location of interest will be a composite signal, and estimates of D H using Eqs.( 36) and (38) will then be estimates of an effective hydraulic diameter, D H,e , that is some function of the different size flow paths that the water traversed.
If the pulse undergoes little modification in shape or duration as it flows through different conduit segments, then the total transmission factor, F T , in a conduit with multiple segments with different values of D H , is given by where F i is the transmission factor from segment i. Furthermore, the total retardation of the thermal peak, τ T , is given by where τ i is the retardation from segment i. τ i is given by where the quantity in square brackets is approximately constant and L i , V i , and D H,i are the length, velocity, and hydraulic diameter, respectively, of segment i.It follows that We can define an effective length (L e ), velocity (V e ), and diameter (D H,e ) such that From this we can see that, provided the quantities in the square bracket are constant, the response of a multi-diameter conduit is the same as that of an equivalent single-diameter conduit with the effective length, velocity, and diameter.We can then consider the relationship of D H,e to D H,i as well as the relationship of L e to L i .There is more than one equivalent conduit that can be defined, depending upon the constraints chosen.We impose the following four constraints, which we deem to be the most physically meaningful: 1.The retardation of the equivalent conduit is equal to that of the multiple segment conduit, 2. Mass (discharge) is conserved along the multiple segment conduit, 3. The equivalent and multiple segment conduits have the same discharge, 4. The flow-through time of the equivalent and multiple segment conduits is the same, Using these constraints, it is possible to solve for the equivalent diameter and length, and where t ft,i is the fluid flow-through time through segment i.
The relationships between the equivalent model parameters and conduit volumes or flow-through times assume that the relationship between hydraulic diameter and cross-sectional area is fixed.An analogous derivation using transmission factor, F , rather than retardation, τ , yields the same relationships, and therefore the equivalent models for damping and retardation are the same.The equivalent diameter is given by the volume-or timeweighted harmonic mean of the hydraulic diameters of the individual segments.Since the harmonic mean accentuates the smaller values in a set, and is always less than the arithmetic mean, one might think that the smaller diameters figure more heavily in the calculation of an equivalent diameter.However, this effect is offset by the weighting by volume, or, equivalently, flow-through time.diameter conduits will have larger volumes and longer flowthrough times.The effect of the weighting is sufficiently strong that, for two conduit segments of equal length, the equivalent diameter is more heavily weighted toward the larger diameter.Since discharge and flow-through time are fixed, the volumes of the multi-segment and equivalent conduits must be the same.Consequently, the length of the equivalent model is equal to this volume divided by the cross-sectional area of the equivalent model conduit.While one might like to hold conduit length fixed between the multi-segment and equivalent models, this is not possible given the constraints (1-4) used above, and we deem these constraints to be more physically meaningful than holding length constant.This is, however, a somewhat arbitrary choice, and other equivalent models could also be derived.As an example of the relationship between multiple segment and equivalent conduit properties, Fig. 4 shows the ratio of D H,e to the average hydraulic diameter, D H,avg , and the ratio of L e to L 1 + L 2 for systems containing two conduit segments with equal length.
Simulations are run in cylindrical coordinates to test if a conduit with two segments with different diameters and a conduit with a constant effective diameter calculated from Eq. ( 48) produce the same transmission and retardation.We run two example cases of a multiple segment conduit, both of which have two segments with equal lengths and different diameters.In one case D H increases by 20 % halfway down the conduit and in the other case by 100 %.Table 2 provides values of model parameters for each case and the transmission and retardation from each simulated thermograph.For the simulations of a conduit with two different D H segments, the output of the first section was used as input into the second section.For the two example cases, there is good agreement between the transmission factors and retardation observed in the multi-segment and equivalent models (Table 2).

Variable flow velocity
During recharge events, discharge variability causes variations in flow velocity, V .To explore the effect of varying velocity on the amount of damping and retardation that occurs, we run additional simulations in cylindrical coordinates where velocity was varied and compared them with constant velocity simulations.For each variable velocity simulation, both V and T w,in are defined by a Gaussian equation of the form of Eq. ( 30).Both curves use the same t peak,in and σ , but the initial velocity, V 0 , and velocity amplitude, V A , (equivalent to T r,0 and R A , respectively, in Eq. 30) are both set to 0.1 m s −1 for all simulations.This simulates a velocity that ranges from 0.1 to 0.2 m s −1 over the duration of the pulse.In these simulations, velocity and input temperature began to change at the same time, and the peak water temperature at the input occurs at the same time as the peak velocity in the conduit.Because these are 1-D simulations of full pipe flow, there are no spatial velocity gradients, even though velocity varies as a function of time.We also run equivalent constant velocity simulations, where the flow velocity for each simulation is set so that the flow-through time is equal to the flow-through time of the peak of the corresponding variable velocity simulation.Five sets of simulations are run to com-pare five different recharge duration, R D , to flow-through time ratios, L/V .
Table 3 provides transmission and retardation data for simulations that consider the effect of variable velocity.For most cases, each set of variable V and constant V simulations produced similar damping.However, as the ratio of recharge duration to flow-through time decreased, the constant V simulations underwent somewhat more damping than variable V simulations.In general, thermal retardation values were similar for the constant and variable V simulations.However, thermal peaks from variable V simulations are characterized by more retardation than the constant V cases when R D L/V and less retardation when R D > L/V .The fraction of time spent at a velocity above or below the average velocity ultimately controlled whether variable V simulations produced less or more retardation, respectively, than the corresponding constant V simulation.The maximum difference between thermal retardation observed in the constant and variable V simulations is approximately 30 %.

Open channel
If water flows along a conduit with a free surface, then a potentially significant amount of heat exchange occurs via radiation through the air.The significance of this exchange is a function of the timescale of the pulse (Covington et al., 2011).To incorporate radiation, we add one more term to the heat advection-dispersion equation: where T s,w is the wet conduit wall surface temperature, h rad is the radiative heat transfer coefficient, A d = P d /W fs is the ratio of dry conduit perimeter (P d ) to the width of the water's free surface (W fs ), A w = P w /W fs is the ratio of wet conduit perimeter (P w ) to W fs , and T s,d is the dry conduit wall surface temperature.h rad is given by where σ SB = 5.67 × 10 −8 W m −2 K −4 is the Stefan-Boltzmann constant.Emissivities of water and rock are close to 1, and temperatures in Eq. ( 51) are in Kelvin.Finally, the dry perimeter boundary condition is where T r,d is the temperature of the dry rock.As before, the wet perimeter boundary condition is given by Eq. ( 8), where T r = T r,w (wet rock temperature), T s = T s,w , and r becomes y.We run three sets of simulations with different choices of recharge duration, R D , with values equal to 1.67 h, 16.7 h, or 6.9 days.For each set, simulations are run in planar coordinates with conduits that are full, mostly full, half full, and mostly empty.A w is held constant for all open channel simulations to see how F and τ vary as a function of A d .All simulations are run with D H = 1 m to further facilitate comparisons.Because ≈ 22 for all of these planar simulations, they accurately model heat exchange in cylindrical or planar conduits and permit simpler planar simulations (Covington et al., 2011).However, we also run a simulation with a full conduit in cylindrical coordinates for comparison to planar simulations for each set.
For the range of A d /A w ratios and R D values considered, there is little difference in the transmission and retardation for each set of simulations with a given recharge duration (Table 4).In general, channels with a free surface undergo slightly more damping than channels that are completely full because there is more rock where heat may be exchanged in the open channel simulations.For the two sets of simulations with longer recharge durations, full planar simulations produce the least retardation, and conduits that are mostly full produce the most retardation.

Thermal recharge shape
Our numerical analysis thus far considers a Gaussian thermal recharge function.This is a rough approximation of the typical shape of thermographs found in natural systems, but natural pulses can display a variety of shapes.To explore the influence of shape on damping and retardation, we run simulations with a variety of other functions that are sometimes used to approximate natural pulses.
of equivalent simulations in cylindrical coordinates.Each set includes a Gaussian function, two types of sine function segments, and a triangular function.Shapes of the recharge thermographs used are shown in Fig. 5.One of the sine-shaped peaks is composed of one period of a sine function from one trough to the next (sine O ) and the other one as half a period between two consecutive zeros of the sine (sine H ). R D for Gaussian functions is 6000 s and 60 000 s, respectively.Sine and triangular functions are defined such that the total area under each curve was equal to the respective areas for the Gaussian functions.In both cases, the sine H curve is the least damped and the triangular thermal recharge is the most damped, although the difference in F between the different recharge functions is small.For the shorter thermal pulse, the Gaussian pulse peaks first, and the triangle, sine O , and sine H peaks occur approximately 4, 4, and 9 % later, respectively, than the Gaussian thermal pulse.For the longer thermal pulse, the triangle pulse peaks approximately 30 % earlier than the Gaussian peak, and the sine O and sine H peaks occur approximately 9 and 17 % later, respectively, than the Gaussian peak.There is less damping and more retardation for thermographs that have a wider peak/trough near the peak/trough, except for the triangle pulse with a R D = 6000 s.However, the triangle pulse is not continuously differentiable, and numerical dispersion likely plays a role.

Luhmann et al. (2012) conducted a multitracer experiment at
Freiheit Spring in southeastern Minnesota by filling a pool next to a sinkhole, heating the pool water, adding tracers, dumping the pool water into the sinkhole, and then monitoring spring breakthrough curves of discharge, temperature, chloride, uranine, delta deuterium, and suspended sediment.The straight-line, horizontal and vertical distances from the sinkhole to the spring are 95 and 19 m, respectively.54 % of the pool's thermal energy was recovered over the first 2 h of the trace, which was lower than either the dye (66 %) or salt recoveries (78 %) over the same time period (Luhmann et al., 2012).The dye recovery was lower than the salt recovery because of degradation, and the lower heat recovery occurred because of the damping of the thermal signal, where some of the heat was transferred into the rock surrounding the flow path.However, the heat from the heated rock was later transferred to subsequent water that flowed along the flow path, since water temperature at the spring remained higher than its background after experiment water no longer reached the spring.The flow path's D H was estimated by reproducing the damped, retarded thermal signal from the trace with heat transport simulations.A much larger diameter was estimated by summing discharge between hydraulic and chemical responses, dividing by flow path distance, and assuming a circular conduit, but the estimate using the thermal pulse was in much better agreement with field observations.We conducted a similar experiment at the same site 3 days later.The pool was filled with approximately 12 600 L of water for the later study.The pool water was heated to 21.5 • C, and 33.02 kg of NaCl were added.Discharge, temperature, electrical conductivity, and suspended sediment data were collected at the spring as the pool water was emptied into the sinkhole.This time, however, the pool was released as two separate pulses.Breakthrough curves are shown in Fig. 6, and all data but suspended sediment time series are provided in Luhmann (2011).Approximately the first half of the 12 600 L of water was released beginning at 16:27 LT on 2 September 2010, and the rest of the pool was emptied into the sinkhole beginning at 16:52 LT.
In general, spring breakthrough curves during this double pulse tracer experiment displayed similar responses to the single pulse tracer experiment 3 days earlier (see Luhmann et al. (2012) for more discussion about the breakthrough curves from the earlier experiment).Discharge at Freiheit Spring increased shortly after each half of the pool was emptied into the sinkhole, suggesting full pipe flow conditions.Furthermore, the initial changes and peaks in suspended sediment occurred before the initial changes and peaks in conductivity.Finally, initial changes and peaks in temperature occurred later than the initial changes and peaks in conductivity because of temperature's non-conservative behavior.
Because these two field-scale experiments were conducted at the same site 3 days apart, all parameters that control F and τ except R D remained nearly constant.There was some rainfall between the two experiments that caused more background variability in spring parameters before the second study, but hydrodynamic conditions were very similar.Background spring discharges before the first and second traces were 26.7 and 26.8 L s −1 , respectively.Additionally, it took 1082 s (Luhmann et al., 2012), 1066 s, and 1103 s between the pool dump (or partial pool dump) and each respective conductivity/chloride increase at the spring for the first trace, the first pulse of the second trace, and the second pulse of the second trace, respectively.Thus, flow-through time was similar for all three pours, and there was little to no variability in D H , L, and V between the two experiments.However, R D was significantly changed from pour one during the first trace (Luhmann et al., 2012) to pours one and two during the second trace.
We did not collect any robust data at the sinkhole during the pours to provide quantitative R D information.However, the time span from the initial increase to the peak in electrical conductivity/chloride at the spring provides a proxy for R D during each pour.This took 625 s during the 30 August 2010 experiment (Luhmann et al., 2012) and 502 and 464 s for the first and second pulses, respectively, of the 2 September 2010 experiment.Because τ is proportional to R 0.5 D , the thermal retardation in planar coordinates of the first or second pulse of the second experiment, τ pl,Ex2 , is given by where τ pl,Ex1 is the thermal retardation in planar coordinates from the first experiment and R D,Ex1 and R D,Ex2 are the recharge durations during single and double pulse experiments, respectively.With τ pl,Ex1 equal to 248 s (Luhmann et al., 2012), the predicted τ pl,Ex2 for the first and second pulses of the second experiment would be 222 s and 214 s, respectively.In reality, τ pl,Ex2 was 224 s and 218 s for the first and second pulses of the double pulse experiment, respectively, providing field evidence that τ ∝ R 0.5 D .Samples were not analyzed for chloride during the double pulse tracer test.Thus, our uncertainty in calculating the transmission factor from either pulse of the double pulse tracer test is larger than our uncertainty from the single pulse tracer test.Furthermore, spring water temperature and electrical conductivity were less stable before the double pulse study because of a recharge event that produced a minimum in conductivity and a maximum in temperature less than 1 day before the beginning of the experiment.Despite these uncertainties, a similar analysis can be performed with transmission using Eq. ( 36), as was done with retardation.The transmission factor in planar coordinates of the first or second pulse of the second experiment, F pl,Ex2 , is given by where F pl,Ex1 is the transmission factor in planar coordinates from the first experiment.With 33 %, respectively.By defining background temperature for each peak as the water temperature before each respective increase in temperature that led to each peak, F pl,Ex2 was 36 % and 34 % for the first and second pulses of the double pulse experiment, respectively.The heat recovery during the first 2 h of the double pulse multitracer experiment (58 %) was higher than the heat recovery over the first 2 h of the single pulse injection (54 %) because of the elevated rock temperatures from earlier water-rock heat exchange.This effect is ultimately responsible for the second pulse producing a higher temperature peak than the first pulse during the double pulse study, even though the second pulse produced a lower conductivity peak with a shorter R D (Fig. 6).The heated rock from the first pulse facilitated the propagation of a higher temperature peak during the second pulse.Thus, while the peak temperature from a later pulse is still useful, deriving flow path information from the peak temperature of a later pulse is more complicated than doing so using peak data from an initial pulse that follows a relatively stable background.
The best simulated fit of the temperature breakthrough curve from the single pulse tracer study occurred with D H = 7 cm using a heat transport simulation in planar coordinates (Luhmann et al., 2012).The average flow-through time between the sinkhole and the spring from the initial increase in discharge to the initial increase in electrical conductivity/chloride at the spring was 1075 s.Given the R D noted above and the values of rock and water physical properties provided in Table 1, then the best D H estimate is 8 cm using τ data from this earlier study and Eq. ( 38).Similarly, the best D H estimate is 5 cm using F data from this experiment and Eq. ( 36).
After these multitracer experiments were conducted, a caver used a track hoe to excavate the sinkhole used for all injections and an abandoned steephead just southwest of Freiheit Spring.Excavation of the sinkhole revealed a relatively flat, weathered bedrock surface with a vertical solution conduit about 20 cm in diameter developed down a vertical joint.The steephead indicates the location of a former spring, which was present long enough for headward erosion to develop the surface feature.The steephead excavation uncovered an underground stream flowing across the back of the steephead toward Freiheit Spring.Although a visual dye trace documented that the steephead flow did emerge at Freiheit Spring, we do not know for sure if water from the multitracer experiments passed through the steephead feature while flowing from the sinkhole to the spring.However, excavation at the steephead revealed a solutionally enlarged bedding plane parting with a height on the order of cms, in agreement with observations at the spring.For a very wide flow channel, D H = 2h, where h is the height of the conduit.Conduit height estimates using either the damping and retardation relationships or the numerical simulations range from 2.5 to 4 cm, in agreement with field observations.9 Discussion

Information content of thermal damping and retardation
Variations in water quantity and quality at karst springs are often used to obtain information about the internal properties of a karst aquifer (e.g., Ashton, 1966;Atkinson, 1977b;Sauter, 1992;Ryan and Meiman, 1996;Birk et al., 2004Birk et al., , 2014;;Luhmann et al., 2011Luhmann et al., , 2012;;Covington et al., 2012).Specifically, Luhmann et al. (2012) showed that combining breakthrough curves of temperature and conservative tracers allows one to constrain values of flow path diameter.This was achieved by adjusting conduit parameters within a numerical transport simulation to obtain best fitting curves for tracer breakthrough.Here, we illustrate an alternative approach that employs the analytical solution for a sinusoidal recharge temperature.This solution provides a good approximation to the damping and retardation of Gaussian temperature pulses simulated over a wide range of conduit properties and recharge conditions.A single fitting parameter, C time , was used to convert between the timescale of the sinusoidal pulse and the timescale of the Gaussian pulse.The primary advantage of this approach is that it is much easier to estimate a hydraulic diameter from analytical equations that relate to damping or retardation than it is to use a numerical model to try to find the best fitting breakthrough curve.Using the technique presented here, one can extract much of the information available in the breakthrough curve using either damping or retardation.
The analytical solution provides explicit relationships for both the transmission (Eq.36) and retardation (Eq.38) of a thermal peak as a function of conduit properties (L and D H ), flow velocity, V , recharge duration, R D , and quantities that are related to the thermal properties of water and rock ( and α r ).The thermal properties of water and rock are relatively constant within a given aquifer, and even do not vary that substantially among near-surface karst aquifers.While an estimate of these parameters is needed to relate damping and retardation to conduit properties, once an estimate is made we typically can treat these as constants for a given site.The conduit length and velocity only occur in Eqs. ( 36) and (38) as a ratio, L/V , which is equal to the flow-through time, t ft .Therefore, we can reduce these two parameters to a single parameter that is also physically meaningful and more easily measured in the field.This leaves three variables, t ft , R D , and D H , that relate to the damping and retardation via two equations.Therefore, if both damping and retardation are measured at a field site, then we have two equations and three unknowns.One might expect that only one of these three unknowns would need to be constrained by additional field data, and then the other two could be calculated from the relations.However, the relations for damping and retardation are not entirely linearly independent, and therefore contain some duplicate information.
A Maclaurin series expansion of the exponential in Eq. ( 36) shows that for low to moderate amounts of damping, the transmission factor, F , scales roughly as Regardless of the extent of damping, the retardation scales as Since t ft and D H enter both relations in the same combination, one of these two variables must be constrained from data in order to solve for the other variables.This conclusion only holds for the low damping regime, but this is also the regime in which damping or retardation could feasibly be measured in the field.These considerations about the independence of the damping and retardation equations are largely theoretical.In realworld cases, both t ft and R D are relatively easy to measure, and it is more likely that both of these will be measured and then used to make separate estimates of D H using both the damping and retardation equations.If these duplicate estimates are substantially different from each other, then it would suggest that some assumptions of the model are being broken or that one or more of the measurements was in error.
Thermal damping and retardation are not affected by the recharge amplitude (R A ) or the thermal conductivity (k w ) or dynamic viscosity of water (µ w ).However, it may be impossible to determine F and τ information if R A is small.Thus, recharge temperatures that are further from background temperatures make it more practical to use water temperature as a tracer to potentially provide flow path information.

Limitations of the model
A fairly large number of simplifying assumptions separate the analytical solution presented above from a natural karst conduit.Therefore, it is worth considering the likely effects of these assumptions, and the extent to which the solution will fail in different settings.Among the assumptions behind the analytical solution are (1) a sinusoidal recharge temperature, (2) a single conduit diameter, (3) no longitudinal dispersion, (4) constant discharge, (5) no hydraulic exchange between the conduit and the matrix, and (6) rock and water thermal properties that are constant throughout the system.
While seasonal temperature variations might be well represented by a sinusoidal solution, most temperature variations at karst springs come in the form of short peaks due to recharge events.However, the numerical simulations presented above demonstrate that, with the help of a fitting parameter, the sinusoidal solution for damping and retardation can be applied to a variety of single-peak functions, including Gaussians, a triangle pulse, and sine peaks.This may not be the case for multipeak functions, particularly if the peaks are more closely spaced than those of the sinusoidal function.In that case, earlier peaks will likely influence the behavior of later peaks.
The analytical solution allows estimation of a single conduit diameter, whereas karst conduits can display a substantial variation in diameter along their length.Therefore, a key question is how this estimated diameter is related to the physical conduit properties.The estimated diameter is the diameter of an equivalent conduit that produces the same thermal damping and retardation.It is possible to derive more than one equivalent model, depending upon the constraints and assumptions applied.However, for a seemingly reasonable set of constraints, the effective diameter is the flow-through time weighted harmonic mean of the hydraulic diameters of the real conduit.To derive this equivalent model, it was assumed that the thermograph timescale does not substantially change as it passes through the system.For the two example simulation sets, the equivalent diameter, so defined, produces the same transmission and retardation as a multi-segment conduit with different diameters (Sect. 7.3.1).This provides some verification that the approach is reasonable, though the approximation is likely to break down for cases where the flow-through time is much longer than the pulse duration.However, this is also the limit in which pulses will be substantially damped and difficult or impossible to observe.
Rather than consisting of a single flow path, natural karst conduits typically contain a network of flow paths of various sizes.Branchwork patterns are quite common, but a variety of network topologies are possible.Physical interpretation of thermal damping and retardation is most straightforward when the system is dominated by a single flow path, such as a sink-rise system.In this case, estimates of conduit diameter apply to the primary conduit.However, a thermal tracing experiment between an injection point and a spring, as conducted by Luhmann et al. (2012), may also allow characterization of a conduit diameter along the flow path between those two points.It is less clear how to interpret natural temperature pulses at a spring fed by a branchwork system, since water arriving at the spring will have flowed via a large number of different paths of different lengths and diameters.In such cases, network properties are likely to play a significant role, and a better understanding of heat transport within networks is required.
The analytical solution also assumes that longitudinal dispersion can be neglected.While karst conduits tend to have high Peclet numbers, and therefore to be advection dominated, dispersion is certain to play a role for increasingly short duration pulses.Therefore, care is needed when applying this solution to short injection pulses, particularly if they propagate a substantial distance.However, the tracer pulses described in Sect.8 are relatively short, and still display the scaling predicted by the theory.In that case, the flow path was also short, which may minimize the influence of dispersion.In addition to longitudinal dispersion, immobile fluid regions such as pools and eddies can substantially influence tracer behavior (Field and Pinsky, 2000).Again, such effects are likely to be largest for short-duration pulses.
The solution assumes constant discharge in time and with distance along the conduit.In Sect.7.3.2,we use simulations to explore the effect of varying discharge in time.We find that discharge variability has a relatively modest effect on damping and retardation, and that the direction of the effect is dependent upon the relative magnitude of the flow-through time and recharge duration.
Our choice of velocity variation is only one of many approximations to natural flow variations, and variability in nature is certainly more complex.While it is difficult to generalize when the constant velocity assumption introduces large errors, the analytical solutions can be used to provide some constraints on the effect of the constant velocity simplification.For example, consider flow through a karst conduit where discharge is not constant, but zero for the first half of the flow-through time and twice the average discharge value for the second half.By comparing to an equivalent conduit with constant discharge and the same water volume in the pulse, the velocity and recharge duration in the non-constant discharge conduit would be 2V and R D /2, respectively.If we determine F , then the term in parenthesis in Eq. ( 36) would be 1/

√
2 times what it would be with constant discharge.Similarly, given τ , the term on the right side of Eq. ( 38) would be 0.5/

√
2 times what it would be with constant discharge.Therefore, the ratio L/(V D H ) in Eqs. ( 36) and ( 38) would be 1/ √ 2 and 0.5/ √ 2 times, respectively, of the actual value.Given an estimate of flow-through time (L/V ), the calculated conduit length would be 1/2 of the true value, and the hydraulic diameter would be underestimated or overestimated by a factor of √ 2 using either the damping or retardation relationships, respectively.Still, it is likely that discharge variability will introduce too much uncertainty and limit the applicability of the damping and retardation analytical solutions in some field scenarios.However, the multitracer studies at Freiheit Spring suggest that the damping and retardation relationships may provide useful results even when there are variations in velocity.
The analytical solution does not account for hydraulic exchange between the conduit and the matrix (or fractures or other conduits).Flow from the conduit to the matrix will affect heat flux in the matrix, and the changed heat flux in the matrix would only have a small, indirect influence on the water temperature in the conduit.In contrast, flow from the matrix to the conduit would have a direct and significant effect on the water temperature in the conduit.However, in this case, the thermal modification of the water is due to mixing of two water sources with different temperatures or dilution of the input rather than damping that occurs due to heat exchange between water and rock.The influence of dilution on the transmission factor can be approximated using a simple mixing model, where the effects of dilution and heat exchange are assumed to be separable.The applicable bounds of this approximation are discussed for linear processes in Covington et al. (2012, Eq. 31), who conclude that the separated treatment of dilution and damping is a good approximation for cases where the signal is not severely damped.The same conclusion applies when heat exchange can be treated as approximately linear, which is also in the regime where damping is not too severe.However, more work is needed to quantify more precisely the conditions under which this separable model breaks down.To account for dilution with a simple mixing model, the peak input temperature is first reduced by the fraction that would be calculated from simple mixing.Then the heat transport model is applied.For example, during the multitracer study in Luhmann et al. (2012), the injected pool water temperature (24.1 • C) produced a peak water temperature at Freiheit Spring of 11.45 • C, above the spring background temperature of 9.08 • C. Without accounting for dilution/mixing, transmission is incorrectly calculated as 16 %.However, chloride concentration from the trace was used to determine the extent of mixing, indicating that the 24.1 • C pool water temperature was reduced to a maximum of 15.19 • C along the flow path due to mixing with water from other sources.By accounting for this mixing, transmission is actually 39 %.Flow from the matrix to the conduit would likely have a small effect on retardation in the conduit that will ultimately depend on the spatial distribution of the matrix input.Further simulations and field experiments could better quantify the effects of dilution/mixing.
Finally, the thermal properties of rock and water are assumed to be constant throughout the aquifer.While the thermal properties of carbonate rocks within karst aquifers can be somewhat variable (Beardsmore and Cull, 2001), uncertainty can be reduced if measured thermal properties for specific formations of interest are available.However, there are still some potential limitations.In particular, many karst conduits contain a substantial layer of sediments on the floor.The heat transfer properties of such sediments are likely to be more variable than that of the solid rock at the field site of interest, and in some cases hyporheic exchange is likely to play an important role.

Considerations for field studies
Determination of the damping and retardation of a thermal peak requires high-resolution data for both temperature and a conservative tracer in order to capture sharp features in the data.In some cases, data output intervals may need to be on the order of seconds to provide sufficient constraints on the timing and magnitude of thermal peaks/troughs.Due to memory or power limitations, data are not often collected at such a high frequency.Consequently, deploying loggers with the capacity to modify data output intervals based on real-time monitoring, or with the capability to transfer data remotely in real time, may be particularly useful.
Monitoring installations in karst frequently have equipment to record water level, electrical conductivity, and tem-perature.In general, water level data have little use in determining retardation, since initial hydrograph perturbations often record arrival of pre-event water.Even in the case of open channel conduits, the discharge pulse, which travels as a kinematic wave, will arrive before the event water.In contrast, spring electrical conductivity perturbations can record event water arrival (e.g., Raeisi et al., 2007), and electrical conductivity interacts more slowly with the rock surrounding a conduit than temperature (Birk et al., 2006;Covington et al., 2012).Thus, in many cases, retardation may be estimated as the time difference between the electrical conductivity and temperature peaks or troughs.
Determining the damping of a thermal peak requires an estimate of recharge temperature, in addition to a thermograph at the spring.In some cases, recharge temperature can be monitored at an upstream monitoring location.If this is not possible, recharge temperature may also be approximated in some special cases, such as during a snowmelt event.Dilution can significantly modify recharge temperatures, and therefore an estimate of dilution is needed for damping calculations, for example by measuring flow at the recharge and discharge points.
While it can be relatively easy to determine thermal retardation using electrical conductivity and temperature data at some monitoring location of interest, interpretation of thermal damping and retardation is most easily accomplished in systems that contain a sinking surface stream.The values of thermal damping and retardation can be estimated during periods of relatively constant discharge between precipitation events.While flow-through time remains relatively constant during these periods, oscillations in surface stream recharge temperature will cause diurnal thermal oscillations at a downstream monitoring location, so long as heat exchange along the conduit is sufficiently ineffective (Luhmann et al., 2011).Measuring discharge at both upstream and downstream monitoring locations allows an estimate of the degree of dilution that occurs along the flow path to facilitate determination of F and constrain potential uncertainty in the measurement of τ and F .Injection of a conservative tracer permits estimates of flow-through time, and thus facilitates calculation of τ when used in conjunction with the travel time of diurnal thermal peaks or troughs from the upstream to the downstream monitoring locations.Measurements of damping and retardation in a sink-rise system are more difficult to obtain during natural recharge events, since temperature and recharge rates may vary independently, and flowthrough time will also vary throughout the event.However, simultaneous monitoring of conductivity and temperature at the recharge and discharge points, particularly if combined with recharge and discharge hydrographs, may still enable measurement of damping and retardation in many settings.
In addition to sink-rise systems, interpretation of damping and retardation may be relatively straightforward during tracer studies with a known recharge input.In this case, the more heavily the system is perturbed, the easier it will be to interpret the results.In general, the ratio of conduit length to the thermal length scale provided in Eq. ( 23) can be used to estimate conditions where it would be possible to perform a thermal tracer study and observe water temperature perturbations at the outlet.This ratio is the thermal process number (Covington et al., 2012), T,sin , and is given by where we use the same relation for ω as in Eq. ( 36).If T,sin 1, then a thermal trace should change water temperature at the outlet, so long as estimates of variables in Eq. ( 57) are appropriate.If T,sin 1, then thermal variations will be completely damped, which still permits estimates of a threshold or maximum conduit diameter (Birk et al., 2014).Regardless of the outcome, thermal tracer studies will generally provide useful results, while either confirming predictions or exposing errors in parameter estimates.
If recharge can be monitored, then R D is given by the full width at half maximum of the recharge thermograph (Eq.33).The actual shape of the pulse will ultimately be a source of uncertainty.When recharge cannot be monitored, a related timescale to the R D is given by the time from the initial change to the peak/trough in a chemograph during a recharge event, as we did in Sect.8.If necessary, the time from the initial change to the peak/trough in a thermograph may be used, although the thermograph will not be as accurate, since the pulse is modified.
Both thermal damping and retardation data can potentially be used to estimate the hydraulic diameter of a karst conduit.However, measurement of retardation, rather than damping, has inherent advantages.There is better agreement in τ between analytical solutions and numerical simulations than there is with F .This suggests that estimates of D H may have less uncertainty when using τ values.Furthermore, it is easier to determine τ in the field than F , since estimates of τ only require temperature and electrical conductivity data at the monitoring location of interest, whereas estimates of F also require information about recharge into the system.Finally, damping requires an accounting of dilute inflow occurring along the flow path.

Conclusions
As water flows through an aquifer, heat exchange occurs between water and rock if they are in thermal disequilibrium.When thermal equilibrium is not attained, the water-rock interaction produces a damped thermal signal in the water that is retarded behind the actual groundwater velocity.Our analytical derivations and numerical simulations demonstrate that the damping (which is quantified using the thermal transmission factor, F ) and retardation (τ ) of thermal peaks in conduits or fractures depend on the flow path's hydraulic diameter (D H ), flow-through time, and the timescale of the temperature variation.Damping and retardation are also dependent on the thermal conductivity, specific heat, and density of rock and the specific heat and density of water.However, these parameters vary relatively little within shallow aquifers.Because of this, the relationships for damping and retardation developed here may be used to estimate the hydraulic diameter of a flow path given estimates of the flowthrough time and the timescale of temperature variations.Our tracer studies at Freiheit Spring provide some evidence for the applicability of these relationships.Additional field work is needed to test the usefulness of these relationships when working with more complex flow paths found in nature.
Simulations with variable D H or velocity, open channels, and sine-or triangular-shaped thermograph shapes produce some variability in F and τ when compared to simulations with constant D H or velocity, full pipe flow, and Gaussian-shaped thermographs.However, variability is generally small, and uncertainty from these conditions should not prevent estimates of D H using F and τ .In general, estimates of D H from natural conduits with variable D H represent a flow-through time weighted harmonic mean of D H .The effect of variable velocity on F and τ relationships is more complex, and additional work is necessary to further understand the effect of the shape and timing of different velocity functions on spring thermographs.Finally, the difference in F and τ between conduits with or without free water surfaces depends on the timescale of temperature variation, but open channels will produce somewhat more damping and retardation than conduits that are water filled.Luhmann et al. (2012) conducted a field tracer experiment that involved temperature, conductivity/chloride, and other parameters.They were able to estimate a flow path's D H using known recharge data, high-resolution output data, and heat transport simulations that reproduced the damped, retarded thermal signal that resulted from the trace.The dependence of F and τ on D H derived here enables a new technique.Specifically, one may estimate the conduit diameter using observations of only the damping and retardation of thermal pulses from natural recharge events or tracer experiments.There is likely more error in D H estimates using this new technique.However, it allows extraction of much of the information carried by the thermal pulses with the ease of employing an analytical solution.Thermal transmission factor in planar coordinates during first pool trace experiment (unitless) F pl,Ex2 Thermal transmission factor in planar coordinates during second pool trace experiment (unitless)

Figure 1 .
Figure 1.Model setup for heat transport simulations involving a (a) conduit or (b) fracture and the surrounding rock.The advectiondispersion equation is solved along the 1-D (a) conduit or (b) fracture.Because of symmetry, conduction in the 3-D rock surrounding the conduit or fracture may be modeled with a simple 2-D rectangle (outlined in thick gray and blue lines).Thus, conduction is modeled in (a) 2-D cylindrical or (b) 2-D planar coordinates.The 1-D conduit or fracture and the 2-D rock are coupled to each other at each respective thick blue line (i.e., the conduit/fracture wall surface).Thick gray limestone boundaries perpendicular to the conduit or fracture are insulated rock boundaries.Thick gray limestone boundaries parallel to the conduit or fracture are sufficiently far from flow path lines to satisfy Eqs.(7) or (9), respectively, and are set to background temperature.
transmission (F planar or F cyl )

Figure 3 .
Figure3.Simulated retardation as a function of theoretical retardation.In general, there is excellent agreement between the analytical solution and numerical simulations.Legend categories are the same as Fig.2, and the graph includes a 1 : 1 line.
Figure 4. (a) D H,e /D H,avg and (b) L e /(L 1 + L 2 ) for different relative increases in D H when L 1 = L 2 .The D H,e for a flow path with two sections of different D H is generally more heavily weighted toward the section with a larger D H , and a larger increase in D H produces a larger D H,e .The L e for a flow path with two sections of different D H is always less than L 1 + L 2 , and a larger increase in D H results in a smaller L e .

Figure 5 .
Figure 5. Different modeled recharge shapes.The sine H curve is widest near the peak and produces less damping and more retardation than the other recharge shapes.Note that the ends of the Gaussian curve are not shown in this figure.
Figure 6.(a) Discharge, (b) suspended sediment, (c) electrical conductivity, and (d) temperature breakthrough curves at Freiheit Spring on 2 September 2010.Water was added to a pool at the edge of a sinkhole, and the water was heated to 21.5 • C as salt was added.The water was emptied into the sinkhole, and breakthrough curves were monitored at Freiheit Spring approximately 95 m away.
planar Thermal transmission factor in planar coordinates (unitless) F T Total thermal transmission factor for multisegment conduit system (unitless) hConduit height (m)h conv Water convection heat transfer coefficient (W m −2 • C −1 ) h rad Radiation heat transfer coefficient (W m −2 • C −1 ) H (1) i Hankel functions of the first kind k r Thermal conductivity of rock (W m −1 • C −1 ) k w Thermal conductivity of water (W m −1 • C −1 )Radial distance from the conduit center into the surrounding rock (m) R Conduit radius (m) R A Recharge amplitude ( • C) R D Recharge duration (s) R D,Ex1 Recharge duration during first pool trace experiment (s) R D,Ex2 Recharge duration during second pool trace experiment (s) Re Reynolds number (unitless) sine H Half period of a sine function between two consecutive zeros (unitless) sine O One period of a sine function from one trough to the next (unitless) t Time (s) t ft Fluid flow-through time through the conduit, L/V (s) t ft,iFluid flow-through time through segment i, L i /V i (s) t peak,in Time of temperature peak at conduit beginning (x = 0) (s) t peak,out Time of temperature peak at conduit end(x = L) (s) T r Rock temperature ( • C) T r,0 Initial rock temperature ( • C) T r,∞Rock temperature at an infinite distance from conduit axis (• C) T r,d Dry rock temperature ( • C) T r,w Wet rock temperature ( • C) T r T r − T r,∞ ( • C) T s Conduit surface temperature ( • C) T s,d Dry conduit surface temperature ( • C or K) T s,w Wet conduit surface temperature ( • C) T w Water temperature ( • C or K) T w,in Water temperature at conduit beginning (x = 0) ( • C) T w,in T w,in − T r,∞ ( • C) T w,out Water temperature at conduit end (x = L) ( • C) Table A1.Continued.T w,out T w,out − T r,∞ ( • C) T w,peak,in Peak/trough water temperature at conduit beginning (x = 0) ( • C) T w,peak,out Peak/trough water temperature at conduit end (x = L) ( • C) V Flow velocity in conduit (m s −1 ) V 0 Initial flow velocity in conduit (m s −1 ) V A Flow velocity amplitude (m s −1 ) V e Equivalent flow velocity in conduit (m s −1 ) V i Flow velocity in conduit of segment i (m s −1 ) V Average or reference flow velocity (m s −1 ) W fs Width of the water free surface (m) x Longitudinal position along conduit (m) y Distance from the conduit center into the surrounding rock (m) α r Thermal diffusivity of rock (m 2 s −1 ) Roughness height (m) Conduction and advection time ratio (unitless) λ T,sin Thermal length scale for sinusoidal temperature variations (m) T,sin Thermal process number (unitless) µ w Dynamic viscosity of water (kg m −1 s −1 ) ρ r Density of rock (kg m −3 ) ρ w Density of water (kg m −3 ) σWidth of thermal Gaussian pulse (s)σ SB Stefan-Boltzmann constant (W m −2 K −4 ) τ Retardation of thermal peak/trough (s) τ iRetardation of thermal peak/trough for segment i (s) τ pl,Ex1Retardation of thermal peak in planar coordinates during first pool trace experiment (s) τ pl,Ex2Retardation of thermal peak in planar coordinates during second pool trace experiment (s) τ planar Retardation of thermal peak/trough in planar coordinates (s) τ T Total retardation of thermal peak/trough for multisegment conduit system (s) ρ w c p,w /(ρ r c p,r ) (unitless)

Table 2 .
Thermal transmission factors and retardation values of variable D H and constant D H,e simulations.
Other parameters different from values in Table1: R D = 6000 s.

Table 3 .
Thermal transmission factors and retardation values of variable V and equivalent constant V simulations.
Other parameters different from values in Table1: L = 5000 m.

Table 4 .
Thermal transmission factors and retardation values of open channel simulations with different A d /A w ratios.
Other parameters different from values in Table1

Table 5 .
Thermal transmission factors and retardation values of different recharge shape simulations.
Other parameters different from values in Table1: L = 5000 m and V = 0.1 m s −1 .

Table A1 .
Notation.A d P d /W fs (unitless) A w P w /W fs (unitless) c p,rSpecific heat capacity of rock (J kg −1 • C −1 ) c p,w Specific heat capacity of water (J kg −1 • C −1 )