Thermal tracer tomography
The most common practice for using heat in subsurface tomography involves stressing the aquifer with a heat signal (Dirac pulse signal, Step signal, etc.) by injecting heated water into the aquifer and moni-toring arrival heat pulses in other locations across the aquifer. One can discern the arrival of heat to a partic-ular location by measuring the temperature with time and establish a temperature breakthrough curve. Breakthrough curves are characterized by their shape (dictated by the heat injection signal), their break-through time (the time when the peak of the breakthrough curve is observed), and the attenuation due to the heat exchanged between heated water and rock’s matrix (de La Bernardie et al., 2018). The transport of heat is linked to the groundwater movement through the preferential paths and the media’s thermal char-acteristics. The temperature evolution curves obtained at different aquifer locations coupled with a mathe-matical model describing the conductive-advective heat transport in subsurface media undergo an inversion process to estimate the aquifer’s hydraulic properties.
The thermal tracer tomography has been the subject of few studies both in porous and fractured media. Doro et al. (Doro et al., 2015) proposed a heat tracer experimental design and its field application to sandy gravel, the alluvial aquifer in Lauwisen Hydrogeological Research Site. They used a dipole forced gra-dient flow and three-level injection system to inject heated water in the middle and ambient water in the above and below sections. They claimed that their results agree well with previous works on delineating the hydraulic properties of the site. Klepikova et al. (2016) also performed a thermal tracer test tomography at an alluvial aquifer in the Hermalle-sous-Argenteau site in Belgium. They also injected a dye tracer simultane-ously with heat and monitored the evolution of temperature and dye concentration in the injection well and closely spaced monitoring wells. They found out that temperature evolution curves were contrary to their expectation of a layered aquifer, explained by unequal lateral and vertical heat transport. They also observed a complex heat plume that is due to lateral and vertical heterogeneity of the aquifer. A temperature-induced water density effect was also noted due to the high temperature of injected water. They finally performed a point-based inversion approach to find the main preferential path by using temperature breakthrough curves. They concluded that thermal heat tracer could be a promising approach for the characterization of subsurface media. In another work, Somogyvari et al. (2016) developed the travel time-based thermal tracer tomography and applied (Somogyvári and Bayer, 2017) this approach to a series of thermal tracer tests per-formed on a shallow alluvial sandy gravel aquifer in the Widen field site in northeast Switzer-land. They per-formed repeated thermal tracer tests with distributed temperature monitoring to obtain multi-source, multi-observation tomographic data. The heated water was injected through a double packer system as a forced hydraulic gradient condition was established. The obtained breakthrough curves and travel time-based to-mography were used to estimate the hydraulic conductivity tomogram. They found that their results agree well with previously reconstructed aquifer hydraulic conductivity tomograms. They concluded that thermal tracer tomography would be a potential approach for resolving the structures and hydraulic conductivity field in heterogeneous aquifers sensitive to the borehole flow. A limitation was nevertheless due to the attenua-tion of the signal that required short distances between the injection borehole and the monitoring boreholes (Somogyvari et al. 2017).
There are not many works reported on using the thermal tracer for subsurface tomography in frac-tured media. Klepikova et al. (2011) presented a tomography approach based on using the fractured bore-holes’ passive temperature profile. They showed a close relationship between the borehole flow profile and the borehole’s passive fluid temperature profile. They used a numerical model to invert the temperature profile to the flow profile. They experimentally validated the proposed approach using temperature data and flow data measured by a heat pulse flow meter. Next, they proposed a three-step inversion approach (Klepikova et al., 2014) to use temperature profiles in the fractured borehole to estimate the transmissivity and connectivity of preferential flow paths in fractured boreholes. In the first step, they determined fracture zones, and then temperature profiles are inverted to flow profiles, and finally, flow tomography is applied to the flow profile obtained in the previous step. They applied the proposed approach to field data obtained in a fractured aquifer in Ploemeur in the northwest of France. They found the calculated fractures transmissivity and their connections correspond well with previous information about the site obtained by geophysical data and other flowmeter tests.
Research motivation and thesis outlines
In the previous section, we show that active FO-DTS is a promising tool that can provide flux meas-urement with great spatial resolution over depth. Compared to the classical hydraulic tomography practices, which rely on a few hydraulic head measurements, A-FO DTS may provide a new interesting source of data that can be useful and sensitive to the transport processes. However, using this new source of data requires an understanding of the advantage and disadvantages compared to hydraulic head data. Chapter 2 describes the advantage of measuring and quantifying the flux data in porous media to reconstruct the aquifer’s het-erogeneity. We use a synthetic numerical model to show that subsurface flux data, which can be measured by active fiber-optic distributed temperature sensors, convey important information about the aquifer’s het-erogeneity. The numerical model is used to simulate the hydraulic tomography tests in which head data and flux data are measured in different locations. Then obtained head data and flux data undergo the inversion process, independently and jointly. The results show that flux data, whether used independently or jointly with head data, can better capture the aquifer heterogeneity than the cases that head data is solely used, especially for a very limited number of data points. As the number of measurements increases, the final inversion result would be independent of the data type.
Chapter 3 presents a real field application of the methodology developed in Chapter 2. The first objective was to test the feasibility of a field experiment designed for hydraulic tomography, including active FO-DTS for measuring fluxes. We thus describe the hydraulic tomography test performed on a shallow sandy gravel aquifer at Saint-Lambert site close to Quebec City. This sandy aquifer has been already used for hy-draulic tomography experiments, and the availabilities of the direct-push methods at INRS makes possible the installation of the FO specifically designed for active DTS measurements. In this hydraulic tomography test, both head and flux data are measured and used to estimate the values and distribution of the hydraulic conductivity. Head data are recorded both in four open boreholes and two boreholes equipped by packers. Two active fiber optic cables buried in the aquifer by the direct push method measure the flux data. It is worth pointing out that to the best of our knowledge, this hydraulic tomography test is the first of its kind and has never been performed. Obtained head and flux data are used within a geostatistical inversion scheme to map the heterogeneity.
It would be interesting to apply the same approach for the fractured media. However, the applica-tion of the direct-push method to deploy the FO cable is currently limited to the shallow unconsolidated aquifer and cannot be used in crystalline rocks. In this case, the only opportunity to measure flux in the fractured media would be to access the media through the boreholes. Maldaner et al. (2019) and Munn et al. (2020) attempted to measure the groundwater flux by using A-FO DTS between the formation and liner. However, the current spatial resolution is too large, for precisely measuring the fracture flow. Furthermore, in the fractured media, the flow varies greatly, and local groundwater flux measured at the borehole scale may not be an appropriate representation of the media’s flow. Thus, in order to counteract these challenges, we suggest measuring the vertical borehole flows using the passive temperature data by extending the work of Klepikova et al. (2011) and also employing the analytical solutions provided by Ramey (1962) and Hassan& Kabir (2002).
Chapter 4 describes how passive heat can allow monitoring the dynamic of subsurface flow in the fractured boreholes. The monitoring process requires real-time measurement of passive temperature data along the borehole attained by fiber-optic distributed temperature sensing. Each recorded temperature pro-file is converted to the flow profile along the borehole, and fracture contributions are determined. This chap-ter shows the added value of using temperature for real-time flow quantification in fractured boreholes, enabling us to capture the dynamic change in the aquifer’s hydraulic state and fractures’ contribution to the total flow.
Finally, Chapter 5 summarizes this work’s key findings and contributions and presents recommen-dations for future researches and development. This thesis comprises five chapters, with three core chapters (Chapters 2, 3, and 4) written as stand-alone manuscripts submitted or in preparation for possible publica-tions in peer-reviewed journals.
Setup of the synthetic test case
To assess the information content of hydraulic head and groundwater flux data for the reconstruc-tion of heterogeneous aquifers, a stationary multi-Gaussian log-hydraulic conductivity field is generated, re-sulting in the field, shown in figure 2-2. The generated aquifer is 550 m in length, 550 m in width and 5 m in depth. The aquifer is discretized into 110 × 110 × 1 in x-, y-, z- directions and corresponding block sizes are 5 m × 5 m × 5 m, respectively. The aquifer is assumed to have one layer (in order to reduce the complexity of analysis of the results and fairly compare the advantages and drawbacks value of each type of data), and the log-hydraulic conductivity field has a multi-Gaussian distribution. The area of interest is chosen in the middle of the aquifer, away from the boundaries, to reduce the boundaries’ effect on the inversion. All boundaries are set to a constant head equal to 350 m.
The correlation length used for generating the Y-field is 75 m and 45 m for x- and y-directions, respectively. The same field (same heterogeneity structure) with different variances of 0.5, 1, 2, and 4 are generated to assess the effect of variance and number of observations on the inversion results. The mean Ymean is -3.57 for all cases, but the Y-fields variance and ranges are different for different experiments.
We use a five-spot setup with a central borehole (P1) and four boreholes (P2, P3, P4 and P5) on the corners of the area of interest (bounded by white dashed lines in figure 2-2). Other monitoring points are also considered between the boreholes, as shown by asterisk symbols. The aquifer is subjected to a series of pumping experiments in each borehole. Two different boundary conditions in the borehole, namely, constant rate and constant head, are considered. When simulating pumping in one borehole, the head and flux values are recorded in other boreholes and monitoring points. The acquired flux and head data are noise-contami-nated before being used to estimate the Y-field.
Hydraulic tomography using head and flux data
Our numerical inversion experiments aim to compare the relative merits of each data type and analyze how borehole boundary conditions (constant head or constant rate), the variance of the Y-field, and the number of observations affect the inversion results. The observational data are generated using the ref-erence hydraulic conductivity field, and a normally distributed error is added to the observations. The stand-ard deviations of these errors are different in all cases. The errors’ standard deviations were chosen to be in a realistic range while ensuring the same initial signal-to-noise ratio of 38 for all cases defined by running the code using Ymean. The resulting measurement errors range from 0.05 to 0.013 (m) for the head and 0.055 to 0.01 (m/day) for the flux.
PCGA with previously mentioned geostatistical parameters are used for the inversion. The trunca-tion order (p-rank) of the prior covariance matrix is chosen as 400 out of 12100. Based on the recommenda-tion by Lee et al. (2016), the truncation order (P), which results in the relative Eigenvalue error below 0.01 would be sufficient to capture most of the covariance matrix structure. The relative Eigenvalue error is de-fined as the ratio of first to (P+1)th Eigenvalue. We have chosen P to be more conservative than having a ratio of 0.01. For Y-field with variance 4, the first Eigenvalue is 1411.47, while the 401st Eigenvalue is 0.047 giving the ratio of 3.25×10-5.
The inversion starts with a constant value of Ymean and continues until the root mean square error between observed and simulated measurement, normalized with the error standard deviation (weighted root mean square error), defined in equation 2:9 reaches a value close to one or after ten iterations. 1 2 WRMSE= √ N∑1 ( – ( )) , N estimate here N is the number of observations and σ is the absolute value of the error’s standard deviation.
Boundary condition at the pumping borehole
Hydraulic tomography is simulated considering two different borehole boundary conditions, con-stant rate (the borehole is being pumped with constant flow rate) and constant head (the head in the bore-hole is kept constant). Note that the external boundary conditions do not change and are kept fixed. The pumping rates used in the model are 2400, 4000, 1750, 5000, and 3800 (m3/day) for P1, P2, P3, P4, and P5, respectively. The equivalent constant head borehole boundary conditions are 324, 340, 300, 329, and 336 (m) for P1, P2, P3, P4, and P5, respectively.
Variance of Y-field and number of observation points
The effect of the Y-field variance is also investigated by considering four different variances (0.5, 1, 2, and 4). Furthermore, a different number of observation points are used to assess their impact on the final inversion results. The observation points are distributed symmetrically in the aquifer. The minimum number of observation points considered is the boreholes (4 observation points) and the maximum are boreholes plus points in the middle shown by asterisk symbols (32 observation points).
The effect of number of observations and variance
Figures 2-7 (a) to (d) show the correlation coefficient (between estimated Y-filed and reference Y-field) versus the number of observations for different type of data and borehole boundary conditions. Blue, red and gray color show the results for head data, flux data and joint inversion of both data, respectively. The data with constant rate borehole boundary condition is marked with solid line while the data with constant head borehole boundary condition is shown by the dashed line.
It can be seen in Figures 2-7 (a) to (d) that as the number of observations increases, the correlation coefficient also increases for all types of data and boundary conditions. For a small number of observations, flux data are superior to head data. The difference between the correlation coefficient of flux and head data is the strongest for a small number of observations, while the difference gradually decreases as the number of observations increases and at a high number of observations, they converge. This is a consequence of the decreasing distance between data points as the number of observations increases, thereby decreasing the radius of averaging. It is worth noting that no added value is observed for joint inversion of both data unless it makes the inversion independent of the type of the borehole boundary condition.
Individual and joint inversion of head and flux data by geostatistical hydraulic tomography The variance of the hydraulic conductivity field affects the final values of the correlation coefficient. The higher the variance, the lower the correlation coefficient (especially for a small number of observations), and also the more challenging it is to reach a WRMSE close to 1. It is worth mentioning that having a higher number of observations, it would be more challenging to converge to the true values.
The effect of truncation order (P) on final inversion results
One of the inversion cases (variance of 4, 32 observation points, joint inversion of head and flux data) was chosen to investigate the effect of truncation order (P) on the final inversion result. Inversions were performed using truncation orders of 25, 50, 100, 200, 400, 800, and 1600. The inversions were per-formed on a server with one Terabyte (1 Tb) memory, 4 processors (Intel Xeon CPU E7-4850 v4 @ 2.10 GHz) and 40 cores in parallel mode. Figure 2-8 shows the Y-field estimated for each P value, the effect on the correlation coefficient, and elapsed time for each geostatistical iteration. For a truncation order of 25, we capture an overly smooth version of the true model with a correlation coefficient of 0.84. By setting the truncation order to 50, 100, and 200, the correlation coefficient increases to 0.88, 0.91, and 0.94, respec-tively. The truncation order of 400 (used in our study) with a correlation coefficient of 0.96 is a turning point beyond which increasing the truncation order does not significantly improve the correlation coefficient. So, the truncation order leads to improvement of Y-field reconstruction up to some points and after this point, it is only the computational time that increases. The computational time increased exponentially for large P values but could be decreased by being less stringent to the final inversion results. The truncation order would be chosen based on the degree of heterogeneity and the computational resource available. It would help perform inversion using a different number of principal components to ensure the proper choice of the number of principal components.
Table of contents :
1.1 General introduction
1.1.1 Subsurface heterogeneity
1.1.2 Tools and techniques for the subsurface characterization
1.2 Inverse problem and hydraulic tomography
1.3 The use of heat for monitoring and characterization of subsurface media
1.3.2 Fiber Optic Distributed Temperature Sensing (FO-DTS)
1.3.3 Borehole flow measurement
1.3.4 Groundwater flow measurement
1.3.5 Thermal tracer tomography
1.4 Research motivation and thesis outlines
Individual and joint inversion of head and flux data by geostatistical hydraulic tomography
2.3.1 Inverse model
2.3.2 Principal Component Geostatistical Approach
2.3.3 Forward model
2.4 Numerical experiments
2.4.1 Setup of the synthetic test case
2.4.2 Hydraulic tomography using head and flux data
188.8.131.52 Boundary condition at the pumping borehole
184.108.40.206 Variance of Y-field and number of observation points
220.127.116.11 Performance Metrics
2.5.1 Inversion of head data
2.5.2 Inversion of flux data
2.5.3 Joint inversion of flux and head data
2.6.1 General findings
2.6.2 The effect of number of observations and variance
2.6.3 The effect of truncation order (P) on final inversion results
2.6.4 Implications for field implementations
Reconstruction of the aquifer heterogeneity using joint inversion of the head and flux data: Application to a shallow granular aquifer
3.2 Experimental setup
3.2.1 State of the art A-FO DTS experiment in shallow aquifers
3.2.2 The Saint-Lambert site
3.2.3 Previous hydraulic tomography experiments
3.2.4 Collection of hydraulic head data
3.2.5 Collection of A-FO DTS data
3.3.1 Inversion of temperature profile to the water flux profile
3.3.2 Hydraulic tomography using PCGA
3.3.3 Groundwater flow modeling
3.4 Field data analysis
3.4.1 Temperature data acquired in FO-17 and FO-18
18.104.22.168 Temperature data quality check and calibration
22.214.171.124 Temperature data to flux profile
3.4.2 Head data in packer-isolated and open boreholes
3.4.3 Setup of the inverse model and hydraulic tomography
3.5 Preliminary results
3.5.1 Hydraulic tomography with pressure data
3.5.2 Hydraulic tomography jointly with flux and pressure data
3.6 Discussion and conclusions
3.6.1 General findings
3.6.2 Field application implications
3.6.3 Perspective works
Modeling borehole flows from Distributed Temperature Sensing Data to monitor groundwater dynamics in fractured media
4.3 Flow model from the temperature profile
4.3.1 State of art
4.3.2 Heat transfer model
4.3.3 Inversion of temperature data and flow profiling
4.3.4 Numerical validation
126.96.36.199 Description of the numerical model
188.8.131.52 Synthetic test with two fractures
4.3.5 Sensitivity of the model to the different parameters
4.4 Hydrogeological setting
4.4.1 The Ploemeur-Guidel field site
4.4.2 Borehole PZ-26
4.5 Data acquisition and processing
4.5.1 Experimental setting
4.5.2 Data processing
184.108.40.206 FO-DTS calibration and processing
220.127.116.11 Heat-pulse flowmeter calibration
4.6 Field application
4.6.1 Spatio-temporal temperature variations
4.6.2 Flow rate calculations
18.104.22.168 Flow Profiling using Distributed Temperature Sensing (DTS) Data
22.214.171.124 Flow Profiling using Heat Pulse Flow meter (HPFM)
4.7.1 Thermal parameters estimation
4.7.2 Added values and drawbacks of using FO DTS for field applications
4.7.3 Inferring groundwater dynamics from temperature data
4.9 Appendix A – Heat Transfer Model
4.10 Supplementary Materials
4.10.2 Calculation Procedure
4.10.3 Automatic flowing zone detection
4.10.4 Linearization of the Ramey’s equation
4.10.5 Inversion of Temperature data and flow profiling
4.10.6 Effect of the different parameters on temperature profile
126.96.36.199 Effect of geothermal gradient on the flowing water temperature
188.8.131.52 Effect of rock thermal conductivity on the flowing water temperature
184.108.40.206 Effect of flow time on flowing water temperature
4.10.7 Heat-Pulse flowmeter calibration curve
5.1 Achieved results
5.2 Future developments