1. Introduction

As an omnipresent parameter, temperature can be seen as a “free tracer” in the hydrosphere. Heat flux mechanisms are indeed showing similarities with solute transport [1]. The acquisition of temperature data has been amazingly useful since decades for numerous and various applications in Earth sciences, to generally infer about complex fluid-flow patterns within geological environments. For example, temperature sensing has been widely used in hydrology for characterizing groundwater-surface water exchange fluxes in streambeds, from temperature sensors [2] , or with optical fiber distributed temperature sensing [3] installed within the hyporheic zone. Shallow groundwater temperature profiles were also investigated to infer recharge conditions, climate and land use changes with use of high-resolution thermistors [4] [5] [6] , level and temperature loggers [7] , thermal analytical modelling [8] , as well as thermal numerical modelling [9]. Depth-temperature logging has extensively been used in hydrogeology, especially for characterizing complex fracture architecture within bedrock aquifers. For that matter, passive measurement with high-resolution thermistors shows great efficiency for mapping hydraulically active fracture distribution and interconnectivity [10] [11] , even by coupling temperature and flow metering surveys together [12] [13]. Advances in groundwater flow characterization in a fractured bedrock aquifer appears to be particularly active using optical fiber DTS, efficiency of which is at its best with an “active temperature sensing” mode, involving the use of a heating device in wellbores [14] [15] [16]. The development of techniques focusing on thermal logging is also remarkable in the field of geological CO₂ storage, where one of the largest concerns is the potential leak of the injected CO₂ from the storage reservoir [17]. Most recent technical developments in this field include permanent monitoring systems installed in injection and monitoring wells, like optical fiber DTS, Fiber Bragg Grating sensors and coaxial cable temperature sensors [18]. Analytical models for inferring CO₂ leakage from temperature logging data have also been recently published [19].

Among other borehole logging techniques, depth-temperature profiling presents the advantages of being a very robust and quick measurement collected with no costly equipment [16]. High resolution temperature logging specifically appears to be extremely sensitive at very low flow [20] , which makes this technique one of the most efficient to sharply characterize very low transmissivity fractures, hydraulic properties distribution in wellbore even at low pumping rate, or to identify very low flow (i.e. passive flows). However, depth-temperature profiles might often appear complex and difficult to decipher for fluid-flow distribution in the wellbore. During pumping, the complexity of depth-temperature signals in wellbore is generally caused by the superimposition of conduction heat flux (i.e. between the aquifer and the borehole wall) with advection heat fluxes (i.e. related to the water inflows distribution). Furthermore, advection heat fluxes relates to the combination of both intensity and temperature of the inflows. This latter characteristic could mislead the direct reading of depth-temperature profiles because temperature shifts measured in the wellbore are then not proportional to the inflow rate. For example, a high inflow rate having quite the same temperature than the one in the water column will produce a small temperature shifts in it. Depth-temperature data are providing sharp and useful information, but the handling of temperature data might be difficult for quantitative interpretation without complementary measurements like flowmetering [12] , or other sophisticated―but time consuming numerical modeling [21].

To deal with these limitations, the objective of this article is to provide a quantitative analytical software that allows the quick interpretation and exploitation of these usually complex distributed depth-temperature profiles in wellbores. Principles implemented in the software calculation rely on a heat budget applied at the scale of the borehole [20] , which is dedicated to the quantification of water fluxes within long screened wellbores. The heat budget is very significantly enhanced in the present work as it is generalized to multidirectional flows in wellbores, allowing to model dynamic depth-temperature data with all possible pumping configurations and fracture distribution with depth. To demonstrate the efficiency of the software, modelling is performed with field data acquired with the best thermistor resolution available (±0.001˚C), and supported by complementary flowmetering and televiewing surveys.

The following abbreviations are used throughout the text for brevity purpose: (STP) static depth-temperature profile of the water column measured in the borehole under “static” condition, i.e. without pumping within the investigated wellbore; (DTP) dynamic depth-temperature profile of the water column measured during the pumping of the investigated wellbore.

2. Materials and Methods

2.1. Study Area

The study area is located in southern Quebec, within the geological region of the St. Lawrence Lowlands (Figure 1). The Ordovician geological units of the St. Lawrence Platform are of sedimentary origin and consist of thick sequences of sandstone of the Potsdam Group, dolomite of the Beekmantown Group, limestone of the Chazy, Black River, and Trenton Groups, Utica shales, and mudstones of the Queenston Group. These geological units are represented in Figure 1 as a simplified version of the detailed mapping by Globensky [22]. The geomorphology of Quebec is marked by glaciation-deglaciation phases, with unconsolidated sediments of glacial and post-glacial origin overlying the fractured bedrock. The complex stratigraphy of the unconsolidated sediment largely controls the hydrogeological context of the underlying fractured bedrock aquifers. In such a glacial geomorphological context, the unconformity between Quaternary unconsolidated sediment and the bedrock is very sharp, and bedrock fracturing generally decreases strongly with depth over the first hundred meters [23].

Borehole loggings were realized during summer 2018 in four long screened wellbores (Figure 1), all drilled in sedimentary bedrocks of the St. Lawrence Plateform, namely in sandstone (wellbore #1), dolomite (wellbores #2 and #4) and limestone (wellbore #3). These wells were accessed in partnership with hydrogeologist consultants who are performing groundwater resource investigation for municipal water supply purpose, as well as access to observation wells of the Ministry of the Environment of Quebec [24]. Positions of the wellbores in Figure 1 are not accurate for confidentiality reasons. All investigated wells have a diameter of 150 mm, are steel-cased within unconsolidated sediments, with casing anchored at least one meter into the bedrock. Boreholes are uncased below the steel tubing. Characteristics of the wellbores are depicted in the synthetic logging figures presented in section 3.1 including: lithology, unconsolidated sediments thickness, steel casing position and length, total borehole depth and depth to the water table (DTWT) measured in static conditions.

2.2. Borehole Logging Materials

“Conventional” logging consisted of measuring water velocities using a spinner flowmeter [25] and performing a camera survey [26] as a support to confirm the nature of active water inflows (Figure 2). Both instruments are trolled into the wellbore using variable size brass centerpieces, and operated with a winch controller connected to a PC [27]. The PC interface allows to know and to adjust the position of the instruments at depth with a resolution of 1 cm, and to control the descent or ascent velocities of the devices into the wellbore in a range of 0.1 to 5.5 m/min. Depth-temperature profiles were measured with a high-resolution thermistor (±0.001˚C) logger equipped with a pressure sensor [28].

The borehole logging surveys were conducted with the same sequence for each investigation, starting with the temperature logging in static conditions, followed by loggings with pumping, including spinner flow metering and transient DTPs measurements. Camera surveys were preferably operated after pumping to ensure the flush of turbid water. Maximum discharge rates were constrained to avoid well dewatering below the base of the steel-casing, allowing temperature, velocity and optical measurements within the whole uncased section of the wellbores. The spinner flowmeter was calibrated for each well under static conditions, with winch down and up speeds varying from 1 to 3 m/min. During pumping tests, the pumps were preferably placed at the top of the wellbore for keeping space for the moving of logging devices through the whole uncased interval of the borehole. Water velocities were logged with the spinner flowmeter trolled downward for upwards flows and vice versa, in order to maximize the flowmeter sensitivity. Flowmeter measurements were performed at a resolution of 5 cm and with a winch down speed of 2 m/min. Raw, noisy signals measured with the flowmeter were smoothed using a moving average of 10 measurements. Flow velocities were converted into flow rates by dividing the measured flow velocities by the section area of the borehole. Total water discharged during pumping was measured with a volumetric counter placed at the hose outlet, and with bucket and chronometer, and compared with the total discharge measured with the flowmeter within the steel-casing. Fluid velocity measurements in the borehole during pumping were taken when steady-state was reached (i.e., with residual drawdown of less than 1 cm/20 min). Pumping tests performed at different discharge rates did not reveal any variation in the

vertical distribution of water inflows into the borehole measurable by the flowmeter. The position of the thermistor at depth during loggings was calculated with the water level measured in static conditions or when water level was stabilized in the wellbore during pumping. When water table was not stabilized during transient temperature loggings, the depth of the thermistor was calculated with data from a pressure logger installed at fixed depth in the wellbore. Sampling frequencies for temperature and hydrostatic pressure were set to one measure per second for all loggings. The thermistor logger was fixed along the cable used for the other devices (spinner flowmeter and camera) allowing it to be centered within the wellbore. Different trolling velocities were tested with the thermistor in stagnant water to ensure that trolling velocities for temperature loggings were not too high, so that the thermistor had sufficient time to equilibrate with its environment when trolled up or down. In stagnant water, trolled velocities up to 4 m/min did not show re-equilibration issues.

2.3. Borehole Heat Budget Calculator

2.3.1. Heat Budget Principles

When a wellbore is pumped, DTP measured in the water column depend on advection and conduction heat fluxes that are in competition. Advection heat fluxes relate to the thermal capacity of groundwater that flows and mixes into the wellbore from hydraulically active fractures. Once pumping is initiated, the borehole wall temperature quickly shifts from the temperature of the geogradient to the temperature of the water flowing in the borehole. Due to this temperature anomaly, conductive heat flux propagates radially from the borehole wall to the surrounding aquifer. A detailed description of principles and heat flux equations involved at the scale of the borehole are given in Meyzonnat et al. [20]. These principles are implemented and generalized in the present development of the present Borehole Heat Budget Calculator (BHB Calculator). The heat budget can be performed in a section of the borehole delimited by couple of temperature measurement with depth, as shown in Figure 3. The general expression of the heat budget provided in this previous work can be rearranged to Equation (1), which provides an explicit incremental expression to model DTPs.

$T_D^z=\frac{q_{in}^{z}C{T}_{in}^z+Q^{z-1}C{T}_D^{z-1}+2dz\frac{\pi \lambda }{\ln \left(\frac{r_e}{r_i}\right)}\left({\bar{T}}_S^{\left(z;z-1\right)}-\frac{T_D^{z-1}}{2}\right)}{\left(q_{in}^z+Q^{z-1}\right)C+\frac{\pi \lambda }{\ln \left(\frac{r_e}{r_i}\right)}dz}$ (1)

With:

・ dz: vertical distance between two temperature measurements in static conditions in the wellbore;

・ ${\bar{T}}_S^{\left(z;z-1\right)}$ (˚C): temperature average of water in the wellbore under static conditions for depth interval (z; z − 1);

・ $T_D^z$ (˚C): the temperature of water in the wellbore at depth under dynamic conditions;

・ $T_{in}^z$ (˚C): the temperature of groundwater inflowing into the wellbore from a fracture at given depth;

・ $q_{in}^z$ (m³/s): the rate of groundwater inflowing into the wellbore from a fracture at given depth;

・ $Q^z$ (m³/s): the water flow rate in the borehole at given depth;

・ C (J∙m⁻³∙˚C⁻¹): specific volumetric thermal capacity of water;

・ r_i (m): the radius of the well;

・ r_e (m): the annular distance of propagation of the temperature anomaly due to pumping into the aquifer;

・ λ (W∙m⁻¹∙˚C⁻¹): the bulk thermal conductivity of the aquifer.

In Equation (1), the bulk aquifer thermal conductivity (λ) is calculated with $\lambda ={\lambda }_s^{\left(1-n\right)}{\lambda }_w^n$ [8] , with aquifer porosity (n), thermal conductivity of bedrock (λ_s Wm⁻¹˚C⁻¹) and thermal conductivity of water (λ_w W∙m⁻¹∙˚C⁻¹). r_e (m) is the time-dependant radius of the heat conduction influence around the borehole. The annular distance of propagation (x = r_e − r_i) of the temperature anomaly (T_D-T_S) in the aquifer occurs as soon as with the beginning of pumping, and is increasing through pumping duration. Given the logarithmic expression for the conduction (Equation (1)), the intensity of conduction fades out with pumping duration (as r_e increases).

In practice, $T_S^z$ are field data measured in static conditions, $q_{in}^z$ , $T_{in}^z$ , and r_e are variables of Equation (1) allowing to model the DTSs measured during pumping. If the flow rates in the water column ( $Q^z$ ) are measured with a flowmeter, water inflows from fractures ( $q_{in}^z$ ) can be calculated, and the heat budget allows to calculate the evolution of the temperatures of water inflows ( $T_{in}^z$ ) during pumping. Such modelling with only temperature logs can provide flow estimation under certain conditions: 1) for extremity flows, the temperature of the first water inflow is known ( $T_{in}^z\approx T_S^z$ ) or when 2) flow rates of the water column ( $Q^z$ ) remain low, so that the dynamic temperature profile provides a unique curvature as advection and conduction remain in competition (see section 2.3.4). For the second case, iterative processing of the heat balance theoretically allows to converge to a single solution for characterizing each inflowing water and to know their respective intensity ( $q_{in}^z$ ) and temperature ( $T_{in}^z$ ) [20].

2.3.2. Generalization of Heat Budget to Any Flow and Pumping Configuration

The modeling of DTPs with the heat budget is incremental and allows to calculate a dynamic temperature at given depth ( $T_D^z$ ), from the previous ( $T_D^{z-1}$ ) using a given direction of water flowing in the wellbore (Equation (1)). Figure 4(a) suggests a configuration with upward flows only. This would be the most common configuration, as in practice it is more convenient to perform pumping tests with a pump placed at the top of the wellbore, so that the whole length of the uncased wellbore is free to be logged for temperature, optical viewing or flowmetering. In this case, heat budget modelling is performed in a context of unidirectional upward flow. However, there are other cases where flow in a wellbore could be oriented downwards, or be even multidirectional within the same well, as represented in Figure 4(b) to Figure 4(e).

In Figure 4(a), flows in the wellbore are all going upwards because the pump is placed at the top of the well. In this context, high-resolution temperature logging should clearly reveal the low inflow from fracture 4, but might not highlight very well the low inflow from fracture 1 because its contribution to the total flow in the upper part of the wellbore is minor. At the opposite, temperature logging in the context of Figure 4(b) (pump placed at the bottom for the same wellbore), should now clearly records the low inflow from fracture 1 but not reveal very well the inflow from fracture 4. In order to reveal both of these low inflows with the same logging, pump could be even placed in an intermediate position (Figure 4(c)), avoiding both of low flows from fractures 1 or 4 to be masked by higher cumulating flows. Some other configuration might also include several pumps placed at different depths of the wellbore (Figure 4(d)). Such a pumping

configuration could be motivated by specific pumping and water treatment operations or because of specific water sampling objectives (i.e. sampling groundwater from fracture 2 only). Several outflows drained out of the investigated wellbore could also be effective through fractures interconnected with a neighbouring pumped well (Figure 4(e)).

The BHB Calculator is implemented in a workbook. In order to generalize the heat budget equation (Equation (1)) towards any possible flow configurations (Figure 4), a series of constrains and logical operations are implemented in the program so that incremental calculations of Equation (1) are automatically applied for a given depth interval with the correct scenario.

2.3.3. BHB Calculator User Guide

The BHB Calculator includes six worksheets. The environment is simplified with editable input data fields colored in blue. All other cells containing calculations are protected and are not visible by users. Results are presented synthetically as spreadsheets and clickable graphics. The first worksheet “Instructions” recalls the information provided in this section as well as the most important concepts regarding the Calculator.

In the “Field data” worksheet, users need to fill out:

・ Properties regarding the wellbore (well radius) and the aquifer (porosity), as well as thermal properties of water and bedrock as defined in section 2.3.1. Typical thermal bedrock properties are suggested at the bottom of the worksheet;

・ STP and DTPs acquired in the field. Each data series may contain up to a maximum of 1200 values. Three different DTP loggings may be input, with specified total discharge rate and time elapsed since the beginning of pumping;

・ Depth-flow profiles measured with doppler or spinner flowmeter devices for each of the pumping tests (three). Each of these flow data series may also contain up to 1200 values.

Control buttons “Clear data” reports to macros that quickly erase long series of data.

The three worksheets “Model” allow to model three DTPs corresponding to three measured profiles filled out in the worksheet “Field data”. The calculator can model any kind of situation regarding positions, intensities of groundwater inflows or outflows (see section 2.3.2). A maximum of 40 different water inflows and outflows can be set for each model. In this worksheet, input data are:

・ Properties of inflows ( $q_{in}^z$ , $T_i^z$ ) and outflows ( $q_{out}^z$ ) and their depth. The depth of each inflow or outflow must be chosen from the available list. The values listed correspond to the depth intervals measured during the STP, which intervals are used for modelling DTPs with the heat budget. Measured DTPs can have different depth intervals than the STP as the model fitting is graphical. Users must enter only one value of inflow or outflow per depth (one per line).

・ Heat conduction settings. These concern depth intervals associated to one value of r_e (see Equation (1)). A maximum of 10 intervals with specific heat conduction can be set. For neglecting the effect of conduction, user may use a high value of r_e (i.e. r_e = 10 m). Conduction intensity is a logarithmic function, so intense conduction is set for r_e very close to r_i, and strongly fades out with increasing r_e.

Warnings to users are programmed to appear and to remain when: 1) the sum of outflows is different from the sum of inflows; 2) r_e input is lower than r_i (provoking a calculation error because of negative parameter for the logarithm); 3) when depth for inflow or outflow is not chosen from one depth interval measured in static conditions; 4) one line contains both inflow and outflow values. Graphics presenting field data (STP, DTPs, and flowmetering data) as well as modeled data (DTPs, flow distribution in the water column) are presented side by side with the spreadsheet containing the model variables, so the effects of variable adjustments on modelled results can directly be visualized.

The fitting procedure (Figure 5) performed by the user consists of adjusting the model variables so that modelled DTPs fits with field data. The fitting procedure must be performed sequentially from the lowest flow (i.e. from extremity inflows) to the highest flows in the water column (i.e. corresponding to the location of the pump). Calculations are incremental as they sequentially follow the flow direction and their increase with depth. Proceeding this way is also an advantage because the temperature of extremity flows can be directly read on field DTP, so it gives solid starting information for the incremental modelling. DTPs usually show unique curvature when both advection and conduction heat flux are in competition [20] , i.e. as long as flow in the water column ( $Q^z$ ) remains low. Under these conditions, DTPs can be modelled with sharp sets of heat budget variables ( $q_{in}^z$ , $T_{in}^z$ , $r_e^{\left[z;z+1\right[}$ ).

If flowmetering data are available, the settings for $q_{in}^z$ , and $q_{out}^z$ directly allow to set the flow distribution ( $Q^z$ ) , so that it corresponds with field data. For

some specific flow conditions or logging operation, the heat budget could self-model DTPs and as such, provide an estimation for both $q_{in}^z$ and $T_{in}^z$ within intervals of the wellbore, or for the whole length of it. See discussion in section 4 and in Meyzonnat et al. [20].

The worksheet “Summary” gives a graphical snapshot of all field and modelled data series. As seen together, results of three measured and modelled profiles might motivate further modeling fitting by adjusting parameters, number and position of inflows, etc. The spreadsheet with modelled data is also provided in this worksheet. Results in the spreadsheet can be selected and copy-pasted. All graphics are editable as users might want to adjust scales for temperatures, depth, or display color and shapes to better represents their data.

3. Results

3.1. Conceptual Modelling Example

The calculator can be used for a wide range of pumping configurations (section 2.3.2). Figure 6 shows the theoretical modelling for three different pumping configurations In this example, sedimentary bedrock is considered, with thermal conductivity of λ_s = 2.5 W∙m⁻¹∙˚C⁻¹ and aquifer porosity n = 0.05. The conceptual example considers a well with five fractures inflows (at depths of 25, 35, 55, 65 and 85 m), a pumping discharge rate of 20 L/min, and a geogradient that is linear as a purpose of simplification. Inflow temperatures are set to the geogradient temperatures, and the intensity of conduction is set at negligible conduction intensity (r_e = 1 m) for flows in the wellbore higher that 2 L/min, and at high intensity (r_e = 0.077 m) for flows lower than that 2 L/min. Values for r_e and for the geogradient slope are arbitrarily and only serves the purpose to provides a simple representation of how the BHB Calculator works. More (complex) cases

are presented in sections 3.2 and 3.3 with in-situ data. The three modelling results in Figure 6 were simply obtained by changing the location of the outflow (the position of the pump) in the BHB Calculator.

Water in the wellbore is flowing upward when pump is placed on top (Figure 6(a)), downwards when the pump is placed at the bottom (Figure 6(c)), and has multidirectional flow directions (with a water direction split) when the pump is placed in the middle (Figure 6(b)). The heat budget Calculator is automatically adapting flows intensities and their directions in the wellbore depending of the position of each inflow and outflow entered by the user. Dynamic temperature profiles are consequently calculated according to each flowing configuration. It is important to notice the curved profiles characterizing lowest flows zones for which conduction is influent (Figure 6). In these theoretical examples, low inflows at the bottom of the wellbore would be better characterized by temperature logging (curved shapes) with pump placed on top, and the low inflow at the top would be better characterized by temperature logging with pump placed on top. The placement of the pump in the intermediate position (between the two higher inflows) would be the most efficient setup for the characterization of the whole inflowing system with only one survey.

3.2. Borehole Logging Results

The STP for borehole #1 (Figure 7) shows a curved shape typical of stagnant water in the wellbore, not disturbed by passive flows in the borehole. The seasonal variation of the temperature at the surface of the soil propagates until 15 m depth and is followed by a reverse geogradient to the bottom of the well. More detailed information about the shape of the geogradient in the context of the Canadian is provided in Gosselin and Mareschal [6] and Meyzonnat et al. [20]. In pumping, the DTPs are very clear and each sudden temperature shift corresponds to the location of a productive fracture. The DTP signal is especially sharp after 55 min since the beginning of the pumping, and highlight water inflows located at 54 m, 44 m, at least two fractures between 35 and 32 m and one weaker but significant inflow at 21 m depth. The spinner flowmeter survey provides the cumulative flow from all active fractures, each flow increment also representing inflow from an active fracture. The spinner flowmeter signal is however affected by noisy signal due to cavity wall effects induced by the variation of the section of the borehole (photos in Figure 7), or when the propeller is disturbed by the (transversal) water inflowing from a fracture. As such, the noisy signal show apparent decreases of the flow before increasing again, which is impossible since in the case of pumping, the water in the wellbore is strictly flowing upwards and flow can only increase with decreasing depth.

As for wellbore #1, the STP for wellbore #2 (Figure 8) shows curved signal characteristic of stagnant water column. However, the seasonal soil temperature variation propagates deeper (i.e. 25 m) than for the previous example, suggesting higher recharges rates. DTPs clearly reveal temperature shifts at 9.1 and 10.9 m indicating major water inflows. These most important water inflows correspond

to the location of small conduits in the bedrock (photos in Figure 8). The DTP measured after 108 min of pumping is most clearly enhancing advection heat fluxes as conduction fades out with pumping duration. On this latter profile (purple color on Figure 8), a small temperature disturbance is distinguishable at 27 m depth facing an impressive cavity. Also, DTPs are not coinciding with the STP for greater depths, indicating the existence a rather low water inflow at 44 m. The measurement with the flowmeter does not highlight any water inflows for depths larger than 11 m. With the latter device, all measurable flows only appear for depth shallower than 11 m. At 27 m, an artefact of apparent drop in flow clearly corresponds to a cavity effect. However, no flow contribution from this cavity is measurable with the flowmeter.

For the wellbore#3 (Figure 9), the STP clearly shows the presence of “induced” flows in the wellbore. At this site, an active municipal well is present at a distance of 200 m from the investigated well and its pumping rates were reported to be as high as 1000 L/min at the time of the field work. From 47 to 70 m, the static profile has a gradient too strong to be natural as for the Quebec context (a gain of nearly 0.6˚C in 20 m instead of being near zero at this depth). The STP shape suggests the presence of water flowing upwards between two (or more) fractures located between 47 and 70 m. With the pumping of the investigated wellbore (pump on top, with flow rates between 25 and 75 L/min), upward flows induced by the neighbouring municipal pumping are further increased. DTPs

confirm upward flows in the wellbore, with rather low flows from 70 to 47 m and a major water inflow at 47 m. The major water inflow at 47 m thus controls the DTPs (vertical profile) until reaching another less productive fracture at 15 m. The flowmeter measurement reports effectively that most of the well’s productivity comes from the cavity present at 47 m depth, as well as the smaller water inflow from cavity at 15 m (photos in Figure 9), but it is inefficient to reveal any flow below 47 m depth.

The STP for wellbore #4 (Figure 10) shows a totally “stepped” signal that instantly indicates the presence of very important flows induced in the wellbore. As for the precedent case, a municipal pumping station is present in the neighbourhood of the investigated wellbore, at a distance of 100 m and water withdrawal reported to be as high as 2000 L/min. The position of hydraulically active fractures at 21, 29, 64, 65, 78 and 88 m of depth are clearly indicated by steps on the STP. Under these conditions, the flowmeter data indicate that the fracture at 88 m (photo at the bottom of the wellbore in Figure 10) is draining out all the water flowing downwards in the wellbore. This bypassed drainage induced by the neighbouring municipal pumping station reaches a total flow of 150 L/min. In pumping conditions (63 L/min, pump placed at the top of the borehole), flowmeter data shows that the water withdrawn at top of the wellbore causes a split of water flow directions at 29.5 m. Above it, upward flow is

discharged by the pump, and below it, downward flows are drained out by the fracture at 88 m. The water split occurs within a major fractured interval between 28.5 and 30 m (top photo in Figure 10). The upper part of the fractured zone drains a mix of water at 9.1˚C, while the lower part of the fractured zone drains a mix of water at near 8.5˚C. The deeper water inflows at 65 and 78 m are cooler and contribute to the cooling of the down flow which is ultimately drained out by the active fracture at the bottom of the borehole.

3.3. Heat Budget Modelling

Heat budget modelling with the BHB Calculator was conducted by following the fitting procedure described in section 2.3.4. If detectable by the spinner flowmeter, water inflow intensities were fixed using these data. Otherwise (for low to very low flow intervals), water inflow intensities were modelled with the heat budget, starting from extremities inflows and subsequently following increasing total flow in the borehole. Such settings were primarily realized for lowest pumping rates having unique DTP allowing to set both inflow temperature and intensity as well as intensity of conduction. The distribution of flow intensities measured with the flowmeter (generally at high pumping rates) were used in heat budgets for lowest pumping rates (changing the total discharge, but keeping the same distribution). For the examples below considering sedimentary bedrock, thermal conductivity of bedrock is set at λ_s = 2.5 W∙m⁻¹∙˚C⁻¹ and aquifer porosity set at n = 0.05.

The thermal modeling for wellbore #1 (Figure 11) was carried out primarily for the smallest flow (22 L/min), 9 min after the beginning of the pumping. Water inflow temperatures at depth were adjusted to the STP ( $T_{in}^z=T_S^z$ ). Then, the conduction was adjusted to fit modeled DTP with those measured in the field (r_e = 0.0754 m). Even for a relatively large flow in the water column above 37 m (10 to 22 L/min), the effect of the conduction towards the geogradient appears not to be negligible. This is due to the high temperature offset due to pumping ( $T_{in}^z-T_D^z$ ) that drives the intensity of conduction. For the most elevated pumping rate (78 L/min) and long pumping time (55 min) it is assumed that conduction intensity had faded out. The DTP that was measured at this pumping rate and duration effectively shows a typical “stepped” signal. For this latter, the conduction was then set to “very low” (r_e = 1 m) and inflow temperatures was set accordingly to match field data. The modeling indicates significant warming of water inflows at 35, 34 and 21 m. and 35 m during the pumping, while the temperatures for the other deeper inflows remained unchanged.

For wellbore #2 (Figure 12), DTPs below 15 m differ from the measured STP. This indicates that there are necessarily water inflows for depths deeper than 15 m, which are however not measurable with the flowmeter. For heat budget modelling, four “deep” fractures (i.e. 44, 40, 35 and 27 m) are placed at the position of slight thermal anomalies visible on DTPs. For shallower depths, the position of three uppermost water inflows (i.e. 12, 13 and 14.9 m) is located at steps

identified with the flowmeter data as well as with DTPs. Heat modeling began with data obtained during the lowest pumping rate (11.7 L/min, duration 39 min). For depths deeper than 15 m, the temperatures of the four water inflows were set equal to the temperature of STP. Then, the intensity of conduction (r_e = 0.0756) and the water inflow rates were adjusted so that the modelled DTPs fit with the in-situ data. For depths shallower than 11 m, it is assumed that higher flows in the water column favours advection, so the heat conduction was set to be weak (r_e = 1 m). With the observed range of flow (i.e. ≥5 L/min for depth ≤ 11 m), even using “high” intensity for conduction has anyway barely no effect on the modelled DTPs. The flow rates of the three upper inflows were adjusted with flowmeter data, and at last, inflow temperatures were adjusted to fit the model. For higher pumping rate and longest pumping duration (39 L/min, duration 108 min), inflow distribution was kept from the first modelling as well as the intensity of conduction. Just by doing so, this second model almost fitted measured

data for depths below 11 m. Adjustment was just made setting a slight warming of the inflow at 27 m. Keeping high intensity for conduction even after relatively long pumping time deepest interval is plausible here because for these depths, flow in the water column remains very weak (i.e. maximum of 1.2 L/min until 15 m depth), so heat advection remains limited. From 15 m depth to the top, the setting of three upper water inflow temperatures was just efficient to fit the model to the data.

For the wellbore #3 (Figure 13), DTPs could be modeled by setting only three water inflows. The shape of the DTP from 10 to 46 m depth is easily interpreted as the consequence of two discrete water inflows at 16 and 46 m. But from 46 to 71 m, the DTP would be either induced by a discrete inflow at 71 m or by a series of very weak distributed inflows. But as there is a temperature shift at 71 m during pumping, the hypothesis of single inflow at 71 m was retained for modelling. Following the same fitting procedure, the first modelling was performed with lowest flow rate (25 L/min) and sufficient pumping time (38 min) so that

DTP is sufficiently distinguishable from the STP. The modelling from 71 to 46 m started by setting the temperature of the bottom inflow ( $T_{in}^z$ ) as the temperature at the bottom of the $T_{in}^z$ log (extremity flow). Then, both conduction intensity (r_e = 0.0755 m) and inflow rate ( $q_{in}^z=0.35$ L/min) were set to fit field DTP. For 46 m to the top of the wellbore, two water inflows were set to correspond to the flowmeter data. Conduction was set to be minor against advection (r_e = 1 m) because of high flow rates (≥20 L/min) for the upper interval. At last, inflow temperatures were adjusted to complete the DTP model fitting. For higher discharge rate (75 L/min, duration 108 min), inflow distribution was kept the same as from the first modelling. Fitting of lower interval required direct setting of the temperature of the extremity inflow, Conduction setting corresponds to slightly less intense value (r_e = 0.0773 m), which makes sense because with increased pumping rate, even the lower inflow is high enough ( $q_{in}^z=1.0$ L/min) to transport heat by advection, participating rather efficiently at the conduction fade out with pumping time. For the upper interval (depth ≤ 46 m), the only changes from modelling at lower pumping rate were to set new inflow temperatures.

For the wellbore #4 (Figure 14), neighbouring municipal pumping induces very high downward flows which are partially compensated by pumping at the top of the wellbore. Since flow rates are high, it is assumed that heat advection dominates against heat conduction in the wellbore, even in the STP. Thus, thermal regime is simple for this example and the main advantage of the modelling is to show the capacity of the calculator to treat multidirectional flows within the same wellbore. Heat modeling consisted of setting inflow and outflow rates from flowmeter measurements, set conduction to be weak (r_e = 1 m) and to adjust inflow temperature to model this “stepped” advection dominated signal. The flow rates set create a split of water flow direction at 29.5 m. Using the heat budget calculator, it becomes obvious that the fractured interval between 28.5 and 30 m actually drains water that does not have the same temperature. Upper water inflows appear warmer and lower water inflows appear colder. Water inflows at

depths of 65 and 88 meters finally bring cold contributions to downward flow at the bottom of the wellbore.

4. Discussions

All temperature profiles surveyed demonstrated the ability of high-resolution DTP to provide sharp information for the localization of water inflows with the flowmeter. The benefit of temperature data quality compared to flowmeter data is visible in DTP for wellbore #1 (depth interval 32 - 35 m and for the inflow at 44 m); wellbore #2 (all signal above 11 m and for the cavity at 27 m), wellbore #3 (inflow at 16 m) and for wellbore #4 (signal around the water split at 29.5 m). Spinner flowmeter is often subject to cavity walls effects (apparent velocity drop within larger irregular section of boreholes), as well as spinner artefacts due to transversal water velocities (rather than vertical) facing water inflow positions. High resolution DTP thus better reports sharp sequences of discrete inflows, even close to each other, where flowmetering signal might be interpreted as “distributed inflows” within one interval. High-resolution DTP surveys also demonstrated very impressive capacity to reveal low inflows. For wellbore #2 (below 11 m depth) and for wellbore #3 (below 47 m), DTPs revealed interesting low-flow patterns which were not detectable with the spinner flowmeter.

Although the BHB Calculator is specifically designed to be used with any possible pump placement, all the DTP in this study were surveyed with the pump placed at the top of the wellbore. This is because the trolling of the thermistor through the wellbore is complicated, or impossible, if cables, tubing and pump(s) are cluttering the section of the borehole. However, depending of fracture distribution in the wellbore, the placement of the pump at depths (i.e. Figure 6(b) and Figure 6(c)) other than at the top of the wellbore (i.e. Figure 6(a)) can certainly improve the pertinence of DTP measurements and their efficiency for quantitative exploitation. For that matter, the use of the BHB Calculator coupled with optical fiber DTS measurements would be better suited for surveys involving pump(s) placed at any intermediate depths in the wellbore. Optical fiber distributed temperature sensing (not used for this study) also certainly provides data with better spatial and temporal resolution compared to thermistors, as it records instantaneous depth-temperature snapshots. But the paradox is that the current temperature resolution of optical fiber might not be still good enough to decipher temperatures ranges in all contexts, especially when collecting temperature data in passive mode and in the context of transitory, near zero geogradient. As a matter of fact, DTPs surveyed in this study with high resolution thermistor clearly showed pertinent temperature patterns (i.e. loggings for wellbores #1 and #2) within a resolution of at least 0.01˚C, which would have certainly not been revealed with optical fiber DTS.

As sharp as they can be, temperature profiles remain difficult to read without the use of the heat budget. Performing quick heat modelling really enhances the comprehension of the temperature logs by clearly visualizing how heat conduction and heat advection can compete depending of flow patterns distributed with the depth. Heat budgeting also removes situations that can lead to confusion, especially when low temperature shifts observed in the wellbore are actually the consequence of high inflow rate, and vice versa. For example, the high inflow at 33 m in wellbore #1 induces a weak temperature shift in the water column because $T_{in}^{33\text{m}}\approx T_D^{33\text{m}}$ . And oppositely, the low inflow at 16 m in wellbore #3 induces a still significant temperature shift in the water column because $T_{in}^{16}\gg T_D^{16}$ .

The fitting procedure of the heat budget with solely temperature data is possible at low flow regime, as long as both heat conduction and advection contributions remain distinguishable in the DTS. For such thermal regime, DTS show curvatures that are best suitable for model fitting. In this situation, the sequential fitting procedure (i.e. starting for extremity flow towards increasing total flow in the wellbore), allows the modelling with constrained values for the parameters ( $T_{in}^z$ , $q_{in}^z$ , r_e). The fitting might be done for the whole length of the borehole with micro-pumping should be operated. Authors demonstrated this latter possibility theoretically with an iterative procedure [20]. When flows into the wellbore become too high, typical stepped temperature profiles indicate that advection dominates. In this case, the heat model cannot be fitted only with temperature data, as one temperature shift in the water column could be modelled by a large range of combinations between temperature or flow rate of each water inflow. High inflow must then be measured with a flowmetering device. There is however a special case for the wellbore #1 survey, since the water inflow temperature corresponds to the temperature measured in static conditions. Temperature of inflows could have been consequently set as values for static temperatures in the heat balance. In doing so, the flows could have been deduced solely by means of the temperature measurement in the wellbore combined with heat modelling.

Another insight provided by heat budget modelling is that temperature of inflows can be calculated. This is where temperature “as a free tracer” might reveal its strongest pertinence among other techniques. The inflow warming or cooling rate during pumping is indicating the orientation of fractures in the aquifer, as well as porosity type (conduit or rather discrete or distributed fracturing). For example, the three upper fractures in wellbore #1 (21, 34 and 35 m) are warming during pumping. This indicates that these fractures are oriented because they necessarily drain water from shallower (warmer) horizons. Such warming of inflows is also suggesting that heat advection is important for these flows within the aquifer, so they could constitute a “conduit” type, more than a dense fractured network. In contrast, deeper inflows in wellbore #1 remain at constant temperature suggesting either drained horizons that have a near-horizontal orientation and/or flows are drained through dense networks of distributed fractures that favor conductive equilibration with the temperature of the bedrock. The same interpretation may be valid for the three upper fractures of wellbore #2 that are warming up during pumping (while those at the bottom do not change). For wellbore #3, all active fractures are cooling down, thus suggesting rather a “conduit” fracturing type, draining deep horizons. For wellbore #4, heat budget modelling makes evident that the fractured interval between 28.5 and 30.0 m does not drain water having the same temperature over the whole interval. Upper water inflow appears warmer and lower water inflow appears colder. This suggests that this less than 2 m thick fractured zone must drain different water reservoirs (a warmer superficial reservoir and a colder one at depth). At last, water inflows for wellbore #4 at depths of 65 and 88 meters finally bring cold contributions that would be associated with deeper reservoirs.

5. Conclusion

Depth-temperature data collected in this study recalled the capacity of the technique to sharply localize water inflows in wellbores. Such profiles collected with a high-resolution thermistor specifically demonstrated tremendous capacity to reveal low flows which remained undetectable with a spinner flowmeter. The Borehole Heat Budget Calculator (BHB Calculator) provided along with this work is dedicated for hydrogeologists who want to enhance the readability and to perform quick quantitative analysis of complex depth-temperature profiles acquired in the context of heterogeneous aquifers. The calculator is easy to use, versatile for heat modeling given any pumping configuration and from any temperature data acquired with optical fiber or thermistor. Optical fiber DTS may provide better spatial and temporal resolution and might be advantageously used for any pump position into the well compared to thermistors. However, thermistors have nowadays much greater temperature resolution the optical fiber, which makes them much more efficient for the characterization of low flows. Explicit modeling of the whole heat system with the BHB Calculator is enhanced in low flow conditions that favor intricate competition between advection and conduction heat fluxes, thus providing unique and explicit curved depth-temperature signals that can be easily deconvoluted. Another insight provided by the use of the BHB Calculator is that the temperature of groundwater inflows can be calculated. This is where temperature “as a free tracer” might reveal its strongest pertinence as it provides information about the origin (i.e. shallower or deeper) of groundwater inflowing into the wellbore. The use the BHB Calculator as an easy tool, combined with the acquisition of temperature profiles using new materials, represents an exciting prospect for the better understanding of flow paths and groundwater origin in complex aquifers.

Acknowledgements

This research was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), through a Collaborative Research and Development grant (CRD) obtained in partnership with TechnoRem Inc. and Envir’Eau-Puits Inc. whose hydrogeological expertise also greatly improved the fulfilment of this research. This work was possible thank to the access to experimental sites that were kindly provided by our partners and by the Direction of Water Policies of the Quebec Ministry of the Environment (MDDELCC).