Travel-time-based thermal tracer tomography

Active thermal tracer testing is a technique to get information about the flow and transport properties of an aquifer. In this paper we propose an innovative methodology using active thermal tracers in a tomographic setup to reconstruct cross-well hydraulic conductivity profiles. This is facilitated by assuming that the propagation of the injected thermal tracer is mainly controlled by advection. To reduce the effects of density and viscosity changes and thermal diffusion, early-time diagnostics are used and specific travel times of the tracer breakthrough curves are extracted. These travel times are inverted with an eikonal solver using the staggered grid method to reduce constraints from the predefined grid geometry and to improve the resolution. Finally, non-reliable pixels are removed from the derived hydraulic conductivity tomograms. The method is applied to successfully reconstruct cross-well profiles as well as a 3-D block of a high-resolution fluvio-aeolian aquifer analog data set. Sensitivity analysis reveals a negligible role of the injection temperature, but more attention has to be drawn to other technical parameters such as the injection rate. This is investigated in more detail through model-based testing using diverse hydraulic and thermal conditions in order to delineate the feasible range of applications for the new tomographic approach.

Main attributes of ideal tracers are their good detectability, their lack of influence on the flow regime, conservativeness and nontoxicity to the environment.Heat is an ideal choice because it is easily detectable by means of traditional temperature sensors, distributed temperature sensors (DTSs) or geophysical techniques (Hermans et al., 2014), and it can be monitored continuously in situ.Typically, background variations are insignificant, and natural heating-cooling cycles have smaller frequencies than the investigated thermal signals.It is also ideal because moderate changes in temperature do not harm the environment, and thus commonly no regulative constraints are imposed.However, due to possible viscosity and buoyancy effects, and their relationship with hydraulic conductivity (K), variation in temperature may modify the flow regime.Ma and Zheng (2010) concluded from numerical simulations that no substantial density effects occur when heating groundwater up by 15 • C.This same critical value is given by Russo and Taddia (2010), based on the recommendations by Schincariol and Schwartz (1990) that buoyancy effects only appear at density differences higher than 0.8 kg m −3 .However, this calculation is only valid if the groundwater temperature is close to 0 • C. By setting a starting temperature of 10 • C (which is more realistic for a shallow aquifer in a temperate climate), this critical density difference is already reached at a heating threshold of 8 • C.This value coincides with that by Ma et al. (2012), who refined their previous findings using field experiments and numerical sensitivity analysis.Essentially, despite several appealing properties, such a tight range for the temperature limits the viability of heat as a tracer.Viscosity and buoyancy effects may render a reliable interpretation of thermal tracer tests impossible.Alternatively, techniques have been developed that can handle broader ranges and are not prone to hydraulic effects of temperature variation.This is the focus of our study.
Our starting point is the fact that for detecting preferential flow paths full analysis of thermal transport behavior may not be necessary.If we focus on characteristic parameters such as travel times or moments of the BTCs, the signal-to-noise ratio may be acceptable for much broader temperature ranges.Travel times of traditional solute tracers are related to the hydraulic properties of aquifers, assuming that the main transport process is advection.This is the case given a sufficient ambient hydraulic or forced gradient during the experiment (Doro et al., 2015;Saar, 2011).One important difference of heat tracer transport over traditional tracers is that diffusion takes place not only in the pore fluid but in the rock matrix as well.So while the tracer front of a solute tracer tends to be sharp, the thermal tracer front appears smoothed.This may make interpretation of BTCs more difficult.
Because thermal diffusion takes place, heat transport is affected not only by the hydraulic properties but by the thermal properties of the aquifer material as well.However, contrasts in thermal parameters are relatively small compared to contrasts in K, which typically spans orders of magnitude (Stauffer et al., 2013).Porosity can also be influential in heat transport, due to the high heat capacity contrast of water and rock components.Yet natural variability in porosity is commonly much smaller than that in K. Therefore variability observed in the transport of a thermal tracer is caused mainly by heterogeneity of K.
In previous studies on thermal tracer testing, diverse setups have been chosen that differ with respect to heating method; injection volumes, rates and temperatures; test duration; and well configurations (Wagner et al., 2014).Mostly hot water is infiltrated in an injection well, and BTCs are recorded in one or more downstream observation well (Ma et al., 2012;Macfarlane et al., 2002;Palmer et al., 1992;Read et al., 2013;Wagner et al., 2014;Wildemeersch et al., 2014).Insight into aquifer heterogeneity is not well constrained by analysis of thermal signals introduced and measured over long screens.To obtain a better definition of the heterogeneity, observations in several wells or at different depth levels need to be compared.Ideally a tomographic setup is chosen, where multiple point injection (sources) and observation points (receivers) are used.By combined inversion of all signals, the spatial variations in K are reconstructed.So far, however, this concept is more established in geophysics and for aquifer characterization in hydraulic tomography, which utilizes pressure signals from depth-dependent pumping tests or multi-level slug tests (Cardiff et al., 2009(Cardiff et al., , 2012;;Illman et al., 2010;Yeh and Zhu, 2007).Klepikova et al. (2014) presented a passive thermal tracer tomography application for characterizing preferential flow paths in fractured media.Their method focused on modeling the fracture network with a sequential method which involves first identifying the location of fault zones on the temperature-depth profiles under ambient flow and pumping conditions.Next, an inversion of the temperature profiles is conducted to obtain borehole flow profiles, and the last step is to estimate the hydraulic properties from these flow profiles.This method provides cross-well connectivities.The work by Doro et al. (2015) is dedicated to the experimental design of cross-well forced gradient thermal tracer tomography.In their approach, a special multi-level injection system is necessary to induce the tracer into a horizontal layer.They also recommend limiting the temperature range to avoid buoyancy effects.Their proposed methodology to interpret the results is to use an inversion scheme developed by Schwede et al. (2014) for this specific experimental setup.This inversion method utilizes the temporal moment of measured BTCs and hydraulic head data together in a joint geostatistical inversion procedure (Illman et al., 2010;Yeh and Zhu, 2007;Zhu et al., 2009).This procedure is computationally demanding, and it assumes a multi-Gaussian distribution of hydraulic properties, which represents a strong restriction in comparison to the true conditions in the field.
In our work, we suggest a travel-time-based inversion procedure, which does not require a priori structural or geostatistical assumptions and is computationally efficient.It is motivated by Vasco and Datta-Gupta (1999), who presented a numerical approach to reconstruct the hydraulic parameters of an aquifer using solute tracer injections in a tomographic setup.As a core element, the transport equation is transformed into an eikonal problem using an asymptotic approach for the tracer transport solution.Their approximation uses the similarity of tracer front propagation to seismic and electromagnetic waves, but with the restriction that the tracer front is abrupt.This approximation can be used for hydraulic signals as well (Vasco et al., 2000), and the travel time of the hydraulic signal can be related to the hydraulic diffusivity of the system.Brauchler et al. (2003) further developed a travel-time-based inversion for Dirac and Heaviside hydraulic sources, using the early-time diagnostics of the signals.To improve spatial resolution, they applied staggered grids (Vesnaver and Böhm, 2000) during inversion.This inversion methodology was applied to several hydraulic laboratory and field experiments (Brauchler et al., 2007(Brauchler et al., , 2011(Brauchler et al., , 2013b;;Hu et al., 2011;Jiménez et al., 2013).Brauchler et al. (2013a) also utilized travel times in a tracer experiment on rock samples on the laboratory scale.Their work revealed that for those samples transport was dominated by the rock matrix, but hydraulic parameters were not estimated.
In this study, we present a new formulation for inversion of spatially distributed hydraulic conductivity using early tracer travel times.It follows the same principles as presented by Brauchler et al. (2003) for hydraulic tomography.Our objective is to obtain a versatile and efficient technique for thermal tracer tomography, which, by focusing on early times, minimizes the role of buoyancy and viscosity effects.In the following section, the new inversion procedure is introduced.It is then applied to a three-dimensional (3-D) high-resolution aquifer analog of the Guarani aquifer in Brazil.We inspect the capability of the new approach to reconstruct 2-D and 3-D sections with heterogeneous K distribution and provide a sensitivity analysis of variable injection rates and temperature ranges.Finally, the findings of exhaustive testing with variable field conditions and technical design parameters are compiled to determine the application window of this new thermal tomography variant.

Travel time inversion
Under high-Péclet-number conditions, when it can be assumed that the thermal transport is dominated by advection, the propagation of an injected thermal plume can be used to gain information about the hydraulic properties of the investigated aquifer.Our goal is to calculate the hydraulic conductivity, K, of the aquifer by inverting the advective thermal tracer breakthrough times.Vasco and Datta-Gupta (1999) showed that the transport equation of a solute tracer can be formulated as an eikonal equation, which is utilized to calculate K.According to this work, a line integral can be written for tracer breakthrough times: ds. (1) Here, t st (x r ) is the breakthrough time of the solute tracer at the receiver (x r ), x s is the source location, v st is the mean tracer velocity, φ is the aquifer porosity and i is the local hydraulic gradient.The line integral relates the tracer breakthrough time to the mean tracer velocity and, thus, to the hydraulic conductivity along the transport trajectory.This equation can be used for a thermal tracer (tt) by including the thermal retardation factor, R: ds. (2) Thermal retardation depends on the porosity of the aquifer, φ; the heat capacity of aquifer matrix, C m ; and the heat capacity of water C w : Changes in these parameters are commonly small compared to changes in K; thus the thermal retardation can be approximated as a constant.For the same reason, φ and the hydraulic gradient, i, are also considered fixed.Values of φ and C can be approximated from prior data, while the hydraulic gradient between observation and injection is measured during the experiment.With these assumptions and the use of standard tomography algorithms, the K distribution can be reconstructed on a pre-defined grid.
In this study, a step function injection temperature signal is used for the active thermal tracer test.In this case the traveling time of the thermal tracer is associated with the propagating thermal front.The tomographic concept requires multiple independent thermal tracer injections at different depths.Temperature BTCs are recorded at multiple observation points, for example at different levels in a downgradient observation well.As common practice for such setups, the number of sources and receivers is one of the important factors that defines the significance and resolution of the results.

Early-time diagnostics
Compared to a conservative solute tracer, heat does not behave ideally.Diffusion is significant in aquifer matrix and pore fluid, while the viscosity and density of the groundwater are variable.Due to the highly diffusive behavior, the emerging thermal front cannot be considered as a sharp transition boundary.In order to obtain accurate results with the inversion, the complications from thermal diffusion need to be mitigated.Both diffusion and mechanical dispersion effects increase with travel time.Mitigation thus can be done by using an earlier characteristic time of the thermal front instead of the (peak of the first derivative) breakthrough time, thus using the fastest component of the heat transport-advection.The earlier characteristic time can then be corrected to the real breakthrough time using a conversion factor, as shown for hydraulic tomography by Brauchler et al. (2003) with a correction for the specific storage coefficient.
The propagation of a thermal front far from the source is described as a one-dimensional (1-D) advection-diffusion problem considering thermal retardation: where R is the thermal retardation factor, T is temperature, D is thermal diffusivity and u is groundwater velocity.
The analytical solution to this problem is (Ogata and Banks, 1961) where T 0 is the initial temperature and erfc is the complementary error function.In this study, we use a step function injection signal as the thermal tracer, and its breakthrough time is associated with the peak of the first derivative of the temperature (Vasco et al., 2000) and can be calculated analytically.During the breakthrough detection, instead of the temperature, the first derivative, T , of the temperature is used as By substituting Eq. ( 5) into Eq.( 6), the peak time can be expressed as Early-time characteristic values can be described proportionally to the peak value: which can be related to the relative peak time (τ α ) as By relating these two expressions, the time of the proportional value can be used to calculate the timing of any value of the signal.Substituting the peak time solution into this expression yields: where f α is the transformation factor that can be used to correct early-time diagnostics back to real breakthrough time.Although Eq. ( 10) has three additional parameters -velocity (u), distance (x) and dispersion coefficient (D) -the function is not sensitive to these values because they are all at higher orders or multiplied with higher orders of velocity.So, by neglecting the terms with higher orders of velocity, they are canceled out.After neglecting the second-order terms of velocity, the expression can be simplified to This equation can be solved analytically for τ α , although infinite numbers of transcendent solutions exist.To have an analytical solution for τ α values between 0 and 1 (times before the peak time), the first branch of the Lambert omega function is applied.The final expression for the transformation factor reads Note that the presented solution is only valid if α and f α are positive.The Lambert omega function is the inverse function of f (W ) = W e W (Weisstein, 2002).Equation ( 13) corresponds to the transformation factor used in hydraulic tomography presented by Brauchler et al. (2003) and Hu et al. (2011).In order to apply the conversion, the temporal scale of the record must be adjusted to the time of the thermal front arrival.In practice, this time is when the first increase on the temperature derivative record can be observed.The application of early-time diagnostics is illustrated in Fig. 1.We are mainly interested in advective transport.However, thermal diffusion may also be significant, smoothening and expanding recorded temperature BTCs, and thus also affecting its derivative.The identification of the peak time through the derivative T is challenging due to the flatness of the curve at the maximum value of the peak.However, using the early-time diagnostics (step 1), only the value of the peak must be known for Eq. ( 8).In step 2, the desired fraction of the peak value (α) and the associated time (τ α t peak ) must be found on the measured T curve.Finally, in step 3, the time is corrected to a calculated peak time using the transformation factor according to Eq. ( 13).In this step, the temperature curve is extrapolated from the fraction time, and by this the effect of diffusion is taken into account.Note that the time zero of the correction is when the thermal front reaches the receiver.This time can practically be chosen when the earliest identifiable temperature change appears at a receiver.
Step 3 allows the travel time to be related to the transport process and to return a real and scaled K value instead of just information about the heterogeneity contrasts.

Staggered grids and null-space energy
To invert the tracer travel times, the SIRT algorithm (simultaneous iterative reconstruction technique) is used to solve the eikonal problem, implemented in GeoTOM3D (Jackson and Tweeton, 1996).The algorithm calculates the transport trajectories between the sources and receivers and solves the line integral of Eq. ( 2) along the trajectories -in a curvebased 1-D coordinate system.To solve the line integral, the solution domain is discretized to a grid.Initially a homogeneous velocity field is defined, and then the velocity values of the cells are updated iteratively to minimize the difference between the inverted and recorded travel times.The algorithm results in mean tracer velocities, and they are transformed into K using the relation of Eq. ( 2), where constant porosity and head gradient are used.In order to provide the uniqueness of the solution, an even-determined problem is needed and thus the number of grid cells should be kept close to the number of measurements (source-receiver combinations).The spatial distribution of the trajectories is never uniform over the domain, the result quality can differ in space and the result can be non-unique (Aster et al., 2011;Menke, 1984).
For discretization, instead of constructing a static regular grid, the staggered grid method (Vesnaver and Böhm, 2000) was used.Solving the problem on a regular grid would highly constrain the freedom of the solution to the geometry of the used grid and the source-receiver locations.By applying the staggered grid method, this constrain can be overcome, with the benefit that the nominal spatial resolution is increased.Otherwise, for a good spatial resolution using one fine grid, a large number of sources and receivers would be required or regularization terms would have to be applied.Staggered grids were successfully employed for hydraulic tomography by Brauchler et al. (2003) and for solute tracer tomography by Brauchler et al. (2013a).In this staggered variant, the problem is solved on different vertically and horizontally shifted versions of a low-resolution regular grid.The inverted results are different for the shifted grids, which are exploited by arithmetically averaging these results to arrive at a final tomogram.The inversion will be stable because of the coarse grids, while the resolution of the averaged tomogram will be as small as the displacements.Although this means that the travel time inversion step will be performed multiple times for one tomogram, it is still computationally affordable due to the marginal computation demand of a single coarse grid resolution.
To characterize the reliability of the results, the null-space energy map is computed.This method has been applied for hydraulic tomography in several studies (Brauchler et al., 2013a, b;Jiménez et al., 2013) and uses the distribution of the inverted transport paths over the inversion grid.The nullspace energy map is calculated from the singular value decomposition (SVD) of the tomographic matrix, which contains the length of each inverted transport path in each grid cell.Values of the null-space energy map are between 0 and 1; thus higher values mean higher uncertainties.Based on the null-space energy map, non-reliable pixels can be deleted from the tomogram.The resulting full inversion procedure, starting with the tracer data and ending with the reliable part of the final K tomogram, is depicted in Fig. 2.

Aquifer analog model
The presented methodology is developed and tested on the Descalvado aquifer analog (Höyng et al., 2014) that is implemented in a finite-element heat transport model (Fig. 3).This analog represents a 3-D high-resolution data set obtained from mapping an outcrop of unconsolidated fluvioaeolian sediments in Brazil.These sediments host parts of the Guarani aquifer system, one of the world's largest groundwater reservoirs.The analog is based on five vertical out-Table 1. Hydraulic conductivity, K; porosity, φ; thermal conductivity, λ; and bulk heat capacity, C; for the nine facies that build up the Descalvado analog.The four zones are introduced for discussion of results and listed here with the major facies components.

Zones
Facies number  crop sections that are recorded during ongoing excavation and interpolated by multi-point geostatistics following the procedure by Comunian et al. (2011).The spatial extent of the analog is 28 m × 7 m × 5.8 m (x, y, z).Hydraulic conductivity, K, and porosity, φ, data were documented on subdecimeter scale, in three parallel and two perpendicular profiles during excavation.Höyng et al. (2014) distinguish nine different hydrofacies (H1-9), which form the primary building blocks and which determine the structural heterogeneity of the characterized volume.In order to ease the interpretation of results, the focus is on major architectural elements, which are the four zones that form the characteristic layers of the formation (Table 1).These can be easily distinguished visually by the dominant color in the selected color scale in Fig. 3: with the blue being top low-conductivity zone, the red central conductive zone, the orange lower-central zone and the yellow bottom zone.In order to use this analog for thermal transport simulations, the original data set (Höyng et al., 2014) is extended with estimated thermal properties (heat capacity, C; thermal conductivity, λ) assigned to the different hydrofacies units ("thermofacies").These properties were calculated based on porosity and available lithological information (Bayer et al., 2015).
The Descalvado aquifer is built up mainly by highly conductive sand and gravel with a layered structure.The average hydraulic conductivity value is approximately K = 10 −4 m s −1 , and the largest difference between two adjacent hydrofacies is three orders of magnitude.Locally, low-K clay intraclasts exist that induce even-higher variations.But, due to sizes of only a few centimeters and a marginal volumetric share, they are negligible for flow and thermal transport simulation.Thermal heterogeneity among the different facies units is controlled by differences in porosity, because the mineral composition does not substantially vary.When clay intraclasts are ignored, thermal conductivity spans from λ = 2.6 to 3.2 W m −1 s −1 , and the volumetric heat capacity ranges between C = 2.4 and 2.6 MJ m −3 K −1 .The global thermal isotropic micro-dispersivity in the forward model is set to about the average grain size, β = 0.1 mm.In the inversion, the mean values were used (β is not used in the inversion because dispersion was neglected).
Flow and transport are simulated as coupled processes, using the software FEFLOW (Diersch, 2014) and the SAMG algebraic multigrid solver (Thum and Stüben, 2012).The analog is embedded into a larger domain with extrapolated homogeneous layers, to minimize lateral boundary effects.The model mesh is generated with the Triangle algorithm (Shewchuk, 1996) and progressively refined towards the analog.Close to wells, the elements are refined to millimeter scale.The total extent of the model is 118 m × 117 m × 15.7 m, consisting in a total of 1 664 626 triangle prism elements.In the center of the model, the resolution of the finite-element mesh is similar to or finer than the resolution of the original aquifer analog data set.
The aquifer is assumed to be confined.In order to simulate initial steady-state conditions with regional groundwater flow in the direction of the long axis, x, constant head boundary conditions are imposed at the perpendicular sides of the model, and no-flow conditions at the other model faces.The constant head values are specified to impose an average hydraulic gradient according to Table 2, but for the inversion the measured cross-well head difference was used.The initial temperature of the model is set to 10 • C.This value is also used as a boundary condition at the sides of the model, which yields isothermal initial conditions.

Experimental setup
We present reconstructions of K fields of 2-D and 3-D analog sections.These sections are called tomograms.Twodimensional profiles represent vertical cross sections between an injection (source) and an observation (receiver) well, while data of three observation wells are utilized for 3-D reconstruction.We specify a base case, which serves as our principal study case, and additionally inspect the perfor-mance of the methodology by varying the experimental design and profile.Note that, independent of the dimensions of the reconstructed sections, the full 3-D analog model was always used to simulate the thermal tracer propagation and resulting travel times, considering buoyancy and viscosity effects.
Focus is set first on 2-D reconstruction.Three profiles in the central plane of the aquifer are selected (Fig. 3).This central plane constitutes a mapped outcrop section with relatively high facies variability.It contains heterogeneous structures of different sizes and contrasts, and it is chosen for being sufficiently far away from the analog boundaries.The location of profile 1 is depicted in Fig. 3. Figure 4a shows the relative locations of an upstream injection well and downstream observation well used for all three 2-D profiles.The distance between the wells is 5 m for an investigated area of 5 m × 6 m.
To examine further the role of aquifer heterogeneity, two additional profiles from the central plane of the analog are investigated.In both cases, the source-receiver geometries are kept the same (Fig. 4a).Profile 2 shows a similar layered structure to profile 1, but with fewer small-scale heterogeneities.The central conductive zone is thicker, providing better connection between the two wells.In profile 3, the central conductive zone is discontinuous, creating a different hydrogeological situation, with weaker connection between the two wells.
In the simulated setup, 6 sources and 6 receivers are employed (Fig. 4a), resulting in a set of 36 source-receiver combinations.The sources are defined as point injections with constant injection rates during the entire simulation time.The used injection temperature signal delineates a Heaviside step function, where the instantaneous change in temperature is arbitrarily set at 0.1 days after the start of simulation, which marks the beginning of the experiment.In order to record BTCs in all observation points even at very small injection rates and temperatures, extremely long simulation times are used (50 days).However, most of the breakthroughs occur during the first five days of the simulation.
The crucial technical design parameters for the experiments are the injection rate, Q, and the injection temperature (or temperature difference, T , in comparison to ambient aquifer conditions).The base values of these two parameters are selected after preliminary field testing (Schweingru-  2 in the sensitivity analysis presented in Sect.4.3. In practice, the source of the injected water can be the investigated aquifer, but note that in this case heating has to be well controlled to keep the injection temperature constant.During a field experiment, the recorded data are always distorted by noise.With the commonly used temperature sensors, this noise is considered very small (Wagner et al., 2014), but still the sensitivity of the temperature sensors is limited.To take this into account when simulating the receiver points, those where the temperature changes are smaller than 0.1 • C are ignored for the inversion.In addition, source-receiver combinations with geometric angles larger than 40 • were not used, following the suggestion of Hu et al. (2011) for hydraulic tomography in layered aquifers.Thus 34 from 36 source-receiver combinations were used in the inversion.
For the 3-D reconstruction, an exemplary case is defined with one injection and three observation wells forming a triangular prism (Fig. 4b) located close to profile 1.The base face is an isosceles triangle, and the observation wells are located along the baseline.The axis of this triangle is at the line where the 2-D profiles are located.The distance between the injection well and the central observation well is 6.5 m, and the length of the triangle base is 3 m.The configuration of the individual wells is the same, resulting in 18 observation points and 108 source-receiver combinations in total.The experiment was simulated using the base values from Table 2, employing the same Heaviside injection signals as in the 2-D cases.

Results and discussion
The following results are structured into four major parts.The first part is the inspection of the inverted tomograms for the three 2-D and one 3-D analog profiles.The second part is the validation of the method using the result of the 3-D reconstruction.The third part is a sensitivity analysis of the inversion procedure with respect to experimental settings such as injection rate and temperature.The fourth part reveals the application window of travel-time-based thermal tomography through rigorous testing with different sections, changing hydraulic conductivity contrasts and varying experimental parameters.

Reconstruction of hydraulic conductivity profiles
The left column of Fig. 5 depicts the analog profiles, and these are contrasted with the inverted ones on the right.For better comparability, the original analogs are upscaled (using the arithmetic mean of the values within a cell) to the same grid as used for the results with 0.125 m × 0.125 m cell size.Figure 5a represents the K distribution of the aquifer analog in profile 1.It is characterized by an overall layered structure, and it shows highest variability with small-scale facies patches in the central part between z = 2 m and z = 4.5 m.Of major interest is the red central conductive zone (hydrofacies H4) at around z = 4 m with non-uniform thickness.In the field, it can cause flow focusing and promote preferential flow.This zone is even more pronounced in profile 2 (Fig. 5c) but not continuous in profile 3, where only laterally high-conductivity wedges can be found.In all profiles, the underlying lower-central zone is dominated by the orange facies H5.With the embedded small-scale layered and crossbedded elements, this zone will give insight into the competence of the inversion procedure to resolve local, decimeterscale structures.
BTCs from 34 source-receiver combinations were used in one tomographic experiment.During staggering, the tomographic inversion is performed on 16 different spatially shifted coarse grids.The uniform cell size of these lowresolution grids is 0.5 m × 0.5 m.In total, 30 iterations are done per inversion, and the inverted velocities are restricted within a range of physically possible tracer velocities.Note that the inversion algorithm allows constraints in velocity to be provided and that, if they are not set appropriately, it can produce outlier pixels close to the sources and receivers, where the flow is focused.Velocity limits (i.e., expected high and low values for K and i) can be calculated using prior information, and the method is not sensitive to small changes in their values.The 16 coarse tomograms are merged together into a fine staggered grid, with a resolution of 0.125 m × 0.125 m.The total computational time for reconstructing one profile was around 10 min on an office PC (Intel ® Core ™ i7-4770 CPU 3.40 GHz).
After calculation of null-space energy maps, a threshold of 85 % is found suitable to constrain the K tomograms.In other words, only pixels with null-space energy of less than 85 % (or vice versa, with a reliability of at least 15 %) are shown in the final reconstructed profile.As illustrated in Fig. 5b, d  and f, this yields fringed edges in the K tomograms and some grayed gaps in the interior.Since the null space denotes local coverage of transport trajectories, there are some regions which are unsatisfactorily accessed.As expected, these are mainly close to the boundaries of the inspected profile and not in the reach of the source-receiver couples.By changing the arbitrary null-space energy threshold, masking of areas of low reliability may be accentuated or mitigated.The most suitable value of the threshold, however, is based on expert knowledge and is set depending on the requirements of the specific case.Experience shows that modifying this value (by 5-10 %) has a minor influence on the visualized structures of major interest, because the null-space energy of the highly conductive zones tends be very small.
The reconstructed profiles in the right column of Fig. 5 shed a first light on the capabilities of thermal tomography.First, we observe that for all profiles the upper zone (in blue) cannot be reconstructed by the inversion.Typically a considerable fraction of it is masked in gray due to the limited contribution to heat transport, which is not surprising due to the low hydraulic conductivity of this zone.In contrast, the tomographic approach identifies the location of the highly conductive upper-central zone (in red) rather well.This zone delineates the fastest travel route between the wells for the heat tracer.Between the upper (blue) and central (red) zones is the strongest contrast in the profiles.This strong contrast shadows the top of the tomograms, because the transport is short-circuited through the high-K zone, with the result that it appears upshifted on the tomogram.When the contrast is smaller, such as in profile 3, this shadow effect is weaker, and it is possible to gain better insight into the low-conductivity zone (Fig. 5e-f).
A striking feature is that the tomographic approach resolves the continuity of the highly conductive upper-central zone in profiles 1 and 2, and it detects the discontinuity in profile 3. Furthermore, the inverted value of hydraulic conductivity of this zone (K = 8 × 10 −4 m s −1 ) is comparable to the original model (K = 1.38 × 10 −3 m s −1 ).For the lowercentral zone, we obtain a similarly good match with an inverted value of 1.6 × 10 −4 m s −1 in comparison to the original value of 2.96 × 10 −4 m s −1 for the dominant hydrofacies H5 (Table 1).This is remarkable, keeping in mind that related travel-time-based techniques of hydraulic tomography have proven to be suited for structural reconstruction, but to a lesser extent for hydraulic parameter estimation (similar match values are found in Brauchler et al., 2007;Cardiff et al., 2013;Jiménez et al., 2013).In many of those studies, parameter values were obtained by ex post calibration with the full forward model (Hu et al., 2011(Hu et al., , 2015;;Jiménez et al., 2013).
The promising findings as depicted in Fig. 5 support the applicability of travel-time-based tracer inversion for thermal tomography, even though thermal diffusion tends to blur ad- vective travel times, which hinders a reliable inversion.However, by taking early arrival times of the recorded BTCs, this effect is minimized.Likewise, when preferential pathways exist, these will be detected by the first thermal breakthrough, which is least influenced by diffusion.As a result, traveltime-based thermal tomography appears especially suited for locating and characterizing high-conductivity zones.
With the 36 source-receiver combinations, exact profile reconstruction is not possible, since the tomograms appear to be smoothed.Fine-scale differences in the form of the highconductivity zone are not reproduced in the tomograms.This is the same for the small facies mosaics that originally occur in the mainly orange lower-central zone.This zone seems mixed with the lower yellow zone, and the hydraulic conductivities of both zones are slightly underestimated.Despite the minor hydraulic contrast between both layers, however, the tomograms indicate locally a facies transition (especially in Fig. 5f).This is not identified in the tomogram of profile 1 (Fig. 5b).Here most small-scale structures exist in the lower-central part above.These cannot be resolved, but they detract from the transport routes of the thermal tracer and thus induce noise in the reconstructions of the lower-central and bottom layer.
Figure 6 shows the reconstruction of the selected 3-D section.The result is presented the same way as the 2-D profiles, using an upscaled version of the original analog for comparison.Three-dimensional staggering is employed, resulting in 64 coarse grids in total.This requires 64 individual inversions and thus a computational time that is drastically longer than in the 2-D cases.With 20 iterations per inversion, the total computational time on the same PC (Intel ® Core ™ i7-4770 CPU 3.40 GHz) was around 1 h for 3-D inversion.The spatial resolution of the coarse grid is 0.5 m × 0.5 m × 0.5 m and of the staggered grid thus is 0.125 m × 0.125 m × 0.125 m.
To assess the reliability of the inverted result, the nullspace energy map is calculated.For the 3-D application a limit of 95 % of reliability is used to accept reconstructed voxels.Lower values would substantially reduce the reconstructed volume, since non-reliable voxels are not presented.Generally, the reliability and thus overall result quality of the 3-D analysis is worse than for the 2-D cases.This is due to the fact that the inverted transport paths cover less of the domain of interest.
Figure 6a depicts the upscaled analog model, sliced in half at the central plane where the injection well is located.The same method of presentation is used for the reconstruction in Fig. 6b.To highlight the differences to the 2-D results, the inverted high-K zone is presented for the whole domain without slicing it in half.The central slice of the 3-D reconstruction is similar to profile 1, because the injection well is located at the same location (Fig. 5a) and the observation wells are located only 1.5 m further away.However, when blanking unreliable voxels in the 3-D visualization, it is difficult to compare the 2-D and 3-D reconstruction in Figs.5b and 6b.At first sight, the reconstructed features of the 3-D and the 2-D inversion are similar.A pixel-to pixel comparison using the central plane of the 3-D reconstruction shows that the difference to the reconstructed values of profile 1 is less than 30 %.This demonstrates that, especially for systems with mainly horizontal structures such as the sedimentary aquifer here, results in 2-D are only slightly improved in a 3-D inversion.Comparing the full profile, the inverted K values are lower than in the 2-D cases but still of the same magnitudes as the original values of the aquifer analog (central conductive zone: 3 × 10 −3 m s −1 inverted to 1.4 × 10 −3 m s −1 original; middle zone: 1 × 10 −4 m s −1 inverted to 3 × 10 −4 m s −1 original).
In Fig. 6b, the central conductive zone of the aquifer is localized mainly at the lateral boundaries close to the wells.Centrally, K values are underestimated and smooth channels appear between injection and observation wells, delineating the suspected main transport paths of the tracer.Similar to the 2-D reconstructions, the central part of these channels is vertically upshifted.The top low-K zone is not reconstructed, but fragments of it appear in the results, marking the location of the contrast boundary on the bottom of this zone.The contrast between the two lower zones can be identified laterally but not centrally -same as in the 2-D profile.
Although neither the 2-D nor the 3-D inversion was capable of reconstructing the top low-K zone, the distribution of the reconstructed transport trajectories can be used to identify these locations.Even though revealing more information about these zones is beyond the scope of this work, it would be an interesting aspect to examine in the future.

Validation
For validation, the reconstructed 3-D K field is implemented in a numerical model with the same settings as used for the forward simulations with the original analog data.Here, homogeneous thermal properties are assumed.In total nine observation wells with six observation points in each are used to validate the inverted result (Fig. 4b).A full tomographic experiment is simulated with six independent warm-water injections using the same configuration as the original simulated experiment.The recorded BTCs are compared with simulations with the aquifer analog data set.The differences in the breakthrough times are used for the validation.
Considering the good reconstruction of the high-K zone, which is most relevant for the thermal transport, we can expect that at most of the observation points the difference would be small.This is exactly what Fig. 7 shows, where the distribution of the differences is presented as a histogram.Most of the values are close to zero, showing a good validation of the result.There are two groups of outliers marked in yellow.The negative outliers are associated with the obser-vations in the top low-K zone where the inversion was not sufficient.Here the predicted heat transport is faster than in the aquifer analog.The second outlier group is related to the underestimated K of the lower-central zone (Fig. 3).The difference in the breakthrough times becomes most significant at observation points that are furthest from the injection well.

Role of injection rate and temperature
The experimental setup may be crucial for the quality of the inversion results.For example, it is well known from related tomographic inversion studies that the feasible resolution depends on arrangement and the numbers of sources and receivers (Cardiff et al., 2013;Paradis et al., 2015).Here we focus on two technical design parameters, which are particularly crucial for thermal tomography when using heated water: the injection temperature and the injection rate.In the following sensitivity analysis, we question whether these need to be carefully tuned or not.Profile 1 is chosen for investigation, depicted again in Fig. 8a and 9a.Note that for forward simulation of travel times the full 3-D analog model is always used.
We first inspect the role of the temperature of the injected water.In all of our models, the ambient groundwater temperature is considered uniform and 10 • C. Viscosity and density effects increase with the temperature difference, T , in comparison to the ambient groundwater.These effects may distort the results of inversion, and thus a maximal difference of T = 8-15 • C has been suggested for thermal tracer testing (Doro et al., 2015;Ma and Zheng, 2010;Russo and Taddia, 2010).This severely constrains the applicability of heat as an active tracer, because it complicates interpretation of BTCs influenced by buoyancy forces.For our tomography, we examine a T from 5 to 80 • C to cover the full range of technical possibilities.The injection rate is kept at Q = 1 L s −1 .
Figure 8 depicts the inverted K tomograms for T = 5, 10, 20, 40 and 80 • C. The results show that the inversion method is not very sensitive to T .The tomograms slightly vary, but they all maintain the major features, and especially the central high-conductivity zone is identified similarly in all variations.Even with an extreme value of T = 80 • C, no distortion appears.This is surprising because buoyancy effects are significant under such conditions.This is attributed to the use of early-time diagnostics, which are mainly controlled by advective transport even if substantial thermal and density gradients prevail in the aquifer.The small differences in the K values can be explained by the changes in viscosity due to the heating.Being able to inject water with high temperature is considered advantageous, because this means that a strong signal is introduced, a high signal-to-noise ratio can be achieved and a greater aquifer volume can be accessed.In practice, of course, maintaining a constant injection temperature at high temperatures can be a technical challenge and requires more sizeable heating devices.
The sensitivity of the injection rate, Q, is investigated in a range of four orders of magnitude, Q = 10 −3 , 10 −2 , 10 −1 , 1 and 10 L s −1 (Fig. 9).The injection temperature is fixed at T = 20 • C. At small injection rates, the heat introduced to the aquifer is small; hence there is no detectable breakthrough at most of the observation points.As shown in Fig. 9b, little insight is obtained with Q = 10 −3 L s −1 , and the quality of the results is poor.Increasing the injection temperature can improve the quality of the result in this case.
By raising the injection rate, the reconstructed continuity of the central conductive zone improves (Fig. 9c-e).For our particular case, this is attributed to the setup.Since the top two observation points are located in the upper lowconductivity zone, this influences the reconstruction of the central high-conductivity zone.
In contrast, at the highest simulated injection rate of Q = 10 L s −1 , the derived tomogram is unsatisfactory close to the injection well (Fig. 9f).This is caused by the highly distorted flow field.Our inversion procedure is based on the assumption that the hydraulic gradient between the two wells is constant.This is not valid anymore, and the relation between inverted mean tracer velocity and hydraulic conductivity is not linear.This effect appears only at very high injection rates, in this case at Q = 10 L s −1 , which exceeds technical possibilities (with an injection temperature of T = 20 • C this would mean 840 kW of thermal power for the experiment).The intensity of the effect of Q settings varies between the different zones.For instance, the lower part of the tomograms in Fig. 9 is not affected.

Application window
The insight gained from variable injection rates and temperatures revealed that the presented tomographic inversion method is robust within a broad range but has limitations.But what exactly are the limits?We tested a broad range of different scenarios to delineate a general application window, where the inversion method can be used to reconstruct the distribution of K in an aquifer.The parameters listed in Table 2 -injection temperature, injection rate and ambient hydraulic gradient -were systematically varied within the given ranges.These ranges were rigorously set, and to reach possible theoretical limits, some scenarios even exceeded the technically feasible range.Additionally, in the three profiles (Fig. 3), the contrasts in the values of K were artificially modified.This was done by expanding or squeezing the original value range for a profile by a factor (range multiplier) between 0.1 and 100.As a result, the original structures of the analog were kept, while the variance was changed.
Each inverted K distribution was compared with the (scaled) analog profile, qualitatively and quantitatively.A first visual test showed whether major structures were reconstructed and the geometries are similar, especially focusing on the conductive zones (Fig. 3).Only acceptable tomograms were kept for the subsequent quantitative analysis.
The quantification is based on an estimated connectivity time between the sources and the receivers.The connectivity time is calculated by converting the K tomogram into a velocity field, using the Darcy equation.With this velocity field, the shortest travel route and time are calculated for all possible source-receiver combinations using the A* pathfinding algorithm (Hart et al., 1968).The root mean square (rms) difference between the connectivity times in the original model and the inverted result is used to quantify result quality relatively to each other and, by this, define an optimal application window for the method.
To condense the results into a normalized parameter space and plot them in a 2-D coordinate system, two dimensionless parameters are selected: the thermal Péclet number (Pe t ) to characterize the hydraulic conditions of the subsurface and the effective injection power to describe the used technical parameters of the experiments.Pe t is calculated separately for the four identified zones of the aquifer: where C w is the heat capacity of the water; q is the Darcy velocity; λ is the thermal conductivity; and d is the length scale, which is here set to unity thickness of the aquifer (d = 1 m).
The used technical parameter effective injection power, P , is defined as where the effective injection rate, Q , represents a normalized rate related to prevailing groundwater flow velocity and calculated for the given length scale, d.Note that Pe t and P are not completely independent; using a higher injection rate can increase the Pe t of a zone.Thus, the defined coordinate system is not orthogonal.
After evaluating approximately 100 different experimental scenarios, resulting in over 350 data points, the application window of the method is identified.In Fig. 10  lines mark strict boundaries between feasible and infeasible regions (where beyond the line no reconstruction is possible), and dashed lines denote an approximate boundary where the result quality of tomograms starts to decrease in the lateral direction (relative decrease in result quality).

continuous
If Pe t is below a critical value, the inversion method is not able to provide any hydraulic information for the investigated zone because the assumption that the heat transport is advective is not valid anymore.In this region, the heat transport is governed by thermal diffusion, and no information on K can be extracted from the heat tracer data.These low-K zones do not build on the resulted tomogram but exist only as high null-space areas.A good example of this is the top lowconductivity zone in Fig. 8b-f, which is not reconstructed properly in any of the presented tomograms.Zones characterized by such low Pe t are typically short-circuited via adjacent conductive zones.The critical Pe t number rises nonlinearly with the increase of P .By raising Pe t with higher injection rate, advection can be promoted in these zones.This provides some information for the tomogram, but the flow field is not short-circuited via an adjacent zone (Fig. 9f), yielding a shadow zone (top low-K zone).
At low P , the amplitude of the tracer breakthrough tends to be too small to be measured in enough observation points to successfully perform the K reconstruction.This strict limit for the application window is due to the assumed 0.1 • C limit for temperature measurement accuracy.It can be overcome by increasing the injection rate or temperature.The result quality gradually declines towards high Pe t and high P .This is caused by the distortion of the flow field from high injection rates (see Fig. 9f).Reconstructions, therefore, may still be acceptable beyond the given dashed boundary.Note that in practice this region is infeasible and hence barely relevant.This is because it corresponds to an injection power of 500 kW-1 MW, and thus this region is also technically infeasible or at least not favorable.

Conclusions
Early arrival times of tracer BTCs are specifically suited for identifying highly conductive zones in heterogeneous aquifers.In our study we formulated a procedure for combined inversion of multiple early arrival times measured during cross-well tracer testing.A tomographic setup with multi-level tracer injection and observation was implemented in a model with a 3-D high-resolution aquifer analog, and we examined the capability of the inversion procedure to reconstruct the heterogeneous distribution of hydraulic conductivity.Heat was selected as a tracer, which offers several advantages in comparison to many solute tracers, but its applicability is traditionally considered limited due to the higher diffusion and coupled thermal-hydraulic processes.
It is demonstrated that the tomographic interpretation of heat tracer signals is well suited for characterization of aquifer heterogeneity.By picking early arrival times, the impact of thermal diffusion, buoyancy and viscosity variation is minimized and, in this way, inversion becomes quasiinsensitive to the temperature range.The presented application window of tested parameters of thermal tracer tomography is wide, and it covers three orders of magnitude for thermal Péclet numbers and five orders of magnitude for injection power.A key principle is that the transport in the aquifer is dominated by advection, and injection of hot water causes minor distortion.This can be controlled, for instance, by establishing a forced gradient between injection and observation point by operating an adjacent pumping well.
The travel-time-based inversion is a fast and computationally efficient procedure, which delivers a tomogram in a few minutes with six sources and receivers.It is revealed that not only structures of mainly highly conductive zones could be reconstructed, but also the values of hydraulic conductivity were closely matched.This is appealing, keeping in mind that the presented eikonal inversion is based on a rough approximation of groundwater flow and transport by a wave equation.Yet when close to strong contrast boundaries, the procedure is not able to reconstruct low-conductivity zones due to short-circuit-shadow effects.To reconstruct these hidden features, a further calibration step or additional information would be required.

Figure 1 .
Figure 1.Three steps of applying early-time diagnostics (ETD) on a thermal breakthrough curve (BTC).(1) Identify the peak T value on the recorded BTC.(2) Find the early-time value to the corresponding fraction of the signal.(3) Extrapolate the early time to the ideal peak time using the transformation factor, f α .

Figure 2 .
Figure 2. Major steps of inversion methodology: (a) conceptual setup of thermal tracer tomography, (b) breakthrough time detection using the early arrival times, (c) tomographic breakthrough time data set, (d) inverted tomograms applying the eikonal solver on different shifted grids, (e) high-resolution tomogram after merging the staggered results together and (f) non-reliable pixels masked after null-space energy calculation.
Figure 3. Vertical cross section through the center of the 3-D Descalvado analog data set showing the distribution of hydraulic conductivity (K).H1-8 represent the hydrofacies units (ignoring clay intraclasts).The location of the three 2-D and one 3-D profile is marked with different colors.

Figure 4 .
Figure 4. (a) Simulated experimental configuration and numerical model boundary conditions.The tomographic setup consists of six sources in the injection well and six receivers in the observation well.(b) Setup of the 3-D experiment with one injection and three observation wells.Additional wells used for validation are marked in gray.

Figure 6
Figure 6.3-D distribution of hydraulic conductivity (K): (a) investigated subdomain of the upscaled aquifer analog and (b) reconstructed tomogram with additional contour lines and unsliced high-K zone.
Figure 7. (a) Histogram plot of absolute differences of breakthrough times between the inverted and the original model (192 samples normalized to the mean of the breakthrough times).Yellow color marks the known outliers, such as observation points in the top low-K zone and the far end of the domain.(b) Scatterplot of observed and simulated breakthrough times.

Figure 10 .
Figure10.The proposed application window of the thermal tracer tomography -related to the injection parameters of the thermal tracer test (effective injection power P ) -and the dominant transport process of the aquifer zone (thermal Péclet number Pe t ).If Pe t is below a critical value, the heat transport is diffusion-dominated, and no hydraulic information can be inverted from the tracer travel times.At low injection power, the temperature change at the observation points is below 0.1 • C and no detection is possible.At very high P the high injection rate distorts the flow field and the results.The application window can help to find the ideal injection parameters based on the prior knowledge about the investigated aquifer.

Table 2 .
Parameterization of experimental setups, with base values and minimum-maximum ranges.