Application of Wavelet De-Noising for Travel-Time Based Hydraulic Tomography

Abstract

Travel-time based hydraulic tomography is a promising method to characterize heterogeneity of porous-fractured aquifers. However, there is inevitable noise in field-scale experimental data and many hydraulic signal travel times, which are derived from the pumping test’s first groundwater level derivative drawdown curves and are strongly influenced by noise. The required data processing is thus quite time consuming and often not accurate enough. Therefore, an effective and accurate de-noising method is required for travel time inversion data processing. In this study, a series of hydraulic tomography experiments were conducted at a porous-fractured aquifer test site in Goettingen, Germany. A numerical model was built according to the site’s field conditions and tested based on diagnostic curve analyses of the field experimental data. Gaussian white noise was then added to the model’s calculated pumping test drawdown data to simulate the real noise in the field. Afterward, different de-noising methods were applied to remove it. This study has proven the superiority of the wavelet de-noising approach compared with several other filters. A wavelet de-noising method with calibrated mother wavelet type, de-noising level, and wavelet level was then determined to obtain the most accurate travel time values. Finally, using this most suitable de-noising method, the experimental hydraulic tomography travel time values were calculated from the de-noised data. The travel time inversion based on this de-noised data has shown results consistent with previous work at the test site.

1. Introduction

Hydraulic tomography has proven to be a reliable method to investigate heterogeneity of hydraulic aquifer parameters at a high spatial resolution level [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]. The application of hydraulic tomography is based on two main steps: hydraulic aquifer testing using a tomographical setup and an inversion of the obtained experimental response data [14,16].

For hydraulic testing, multiple investigation application methods such as pumping tests and slug tests may be employed in a setup comparable to computerized cross-section scanning (CAT) or seismic tomography (ST) [17,18]. The key to obtain optimal subsurface characterization results is then the selection of suitable inversion methods to handle the generally huge amount of obtained experimental data.

Hydraulic tomography inversion may be performed using a numerical groundwater flow model and parameter estimation [5,11], providing spatial distribution of hydraulic conductivity (K) and specific storage (Ss) values. This inversion method can provide remarkably accurate inversion results. However, it may be too time-consuming [19,20].

Another hydraulic tomography inversion method is based on the inversion of hydraulic signal travel times, which can provide two- or three-dimensional spatial distributions of hydraulic diffusivity (D) [14,16,21,22,23,24]. This approach is based on the solution of a line integral equation, which relates the travel time of a hydraulic signal between an excitation point and an observation point with hydraulic diffusivity between the points [14,25]. In this way, the groundwater dynamics equation is transformed into a one-dimensional linear integral, which allows users to quickly handle a huge amount of tomographical data on a personal computer. In addition, since the signal travel time may be relatively short, long-time hydraulic testing is often not required, which makes field work less time-consuming.

However, the signal travel time values, which can be very small (i.e., less than 1 s), are strongly sensitive to noise in the experimental data. The noise in experimental data can have a complex structure. A possible source for this may be the inherent noise generated by electronic components, the noise caused by the design error of the circuit, the defects in the installation process, or the interference noise (e.g., spatially radiated interference noise) [26]. Therefore, the de-noising method plays a crucial role when travel-time based inversion is applied. Noise is a general term for unwanted (and unknown) modifications that signal may suffer during capture, storage, transmission, processing, or conversion [27]. It is assumed here that the main source of noise in the experimental data is the inherent noise generated by electronic components. This kind of noise is also known as thermal noise. It is generated by the agitation of charge carriers (usually electrons) inside an electrical conductor at equilibrium, which happens regardless of any applied voltage [28]. It is approximately white in an ideal resistor, which means its probability distribution follows a normal function [29]. Therefore, it can be regarded as Gaussian white noise.

Noise reduction is the recovery of the original signal from the noise-corrupted one and is a common goal in the design of signal processing systems. Various noise reduction techniques were used in multiple application fields. For instance, Illman et al. [30] and Liu et al. [12] used a polynomial fitting technique to fit the drawdown data obtained from sandbox experiments. Wavelets, as a relatively new signal processing tools, have found successful applications in a variety of signal processing problems [31,32], including de-noising. Xie et al. [33] used a wavelet transform to de-noise speckled Synthetic Aperture Radar (SAR) images and achieved better results compared to the images de-noised by a refined Lee filter. Quiroga et al. [34] applied a wavelet de-noising tool to denoise simulated data and achieved a significantly better reconstruction of event-related responses (ERPs) in comparison with a reconstruction based on a conventional Wiener filtering. Xiang et al. [35] applied a wavelet de-noising technique to process the data from sandbox experiments for a SimSLE analysis to map aquifer heterogeneity. However, there are not many previous studies which have applied wavelet de-noising to analyze data from hydrogeological hydraulic field experiments [14,36,37].

Therefore, the purpose of this study is to compare the performance of different de-noising methods to find the most appropriate de-noising method to process data from groundwater hydraulic experiments. A numerical groundwater flow model was built according to the real field situation at our fractured-porous aquifer test site to simulate ideal drawdown data from tomographical pumping tests. Then, Gaussian white noise was added to the ideal data to simulate real noise flawed measurements. Different de-noising methods were then applied to correct the noisy data and hydraulic signal travel time values were calculated with the de-noised data and compared to the ‘true’ travel time values derived from the ideal data.

In order to prove that the numerical model provides similar aquifer conditions and characteristics compared with the real test site, a diagnostic plot analysis method was applied to drawdown data in this study. A diagnostic plot is a simultaneous plot of drawdown data and their logarithmic derivative as a function of time in log-log scale [38,39]. It is demonstrated that the diagnostic plot approach allows for a reliable identification of the characteristics of the aquifer [39], and for a validation of the numerical aquifer model.

2. Methodology

2.1. Hydraulic Travel Time Inversion

The main methodology for travel-time based inversion is based on the findings by Vasco and Datta-Gupta and Brauchler et al. [14,25]. Through their work a relationship between the travel time of a hydraulic signal and the diffusivity within the investigated domain can be established as Equation (1):

$$\sqrt{t_{peak}}=\frac{1}{\sqrt{c}}{{\displaystyle \int }}_{x_1}^{x_2}\frac{ds}{\sqrt{D\left(s\right)}},$$

(1)

where $t_{peak}$ is the travel time of the pressure signal from the source $x_1$ to the receiver $x_2$; $s$ is the propagation path of the signal; $D\left(s\right)$ is the diffusivity along the path; and c is a dimension dependent coefficient. For a two-dimensional (2D) situation it is c = 4; for three-dimensional (3D) it is c = 6.

Figure 1 shows the area Ω between the pumping well and observation well. The area is divided into 15 cells, each denominated discretely as ${\Omega }_j$, where j = 1, …, 15. $D_j$ represents the average value of hydraulic diffusivity within the jth cell.

Equation (1) can be rewritten in discrete form, as [40]:

$$\sqrt{c{t}_i}={\displaystyle \sum }_{j=1}^J\sqrt{\frac{1}{D_j}}{s}_{ij}, i=1, \dots , I, j=1, \dots , J,$$

(2)

where I is the total number of signals, J is the total number of the cells within the investigation are, anc $s_{ij}$ is the path length of the ith signal through the jth cell.

Equation (2) is reformulated into a matrix form as:

$$T=SB$$

$$T=\sqrt{c{t}_i}, S=s_{ij}, B=\sqrt{\frac{1}{D_j}}, i=1, \dots , I, j=1, \dots , J.$$

(3)

2.2. SIRT Algorithm

Simultaneous iterative reconstruction technique (SIRT) is applied for solving Equation (3) to invert hydraulic diffusion coefficients (D) based on travel time data [40,41].

After setting the current $B_j^{est}$ equal to the initial value $B_j^{init}$ for j = 1, …, J, the following are the steps of SIRT to update the $B_j^{est}$ value.

Step 1: Perform calculation with Equation (3), providing a predicted $T_i^{pre}$:

$$T_i^{pre}={{\displaystyle \sum }}_{j=1}^J{S}_{ij}{B}_j^{est}, i=1, \dots , I.$$

(4)

Step 2: Find the correction for each cell by examining the signals passing through that cell and averaging the corrections suggested by each signal. This operation for the jth cell is defined by:

$$\Delta B_j=\frac{1}{W_j}{{\displaystyle \sum }}_{i=1}^I{\Delta }^i{B}_j=\frac{1}{W_j}{{\displaystyle \sum }}_{i=1}^I{S}_{ij}\frac{T_i^{obs}-{\sum }_{j=1}^J{S}_{ij}{B}_j^{est}}{{\sum }_{j=1}^J{S}_{ij}^2}, j=1, \dots , J.$$

(5)

The weight $W_j$ is the number of signals intersecting the jth cell in order to obtain an average correction $\Delta B_j$.

Step 3: Determine the new $B_j^{\left(new\right)est}$ with the average correction $\Delta B_j$:

$$B_j^{\left(new\right)est}=B_j^{est}+\Delta B_j, j=1,\dots ,J.$$

(6)

Then go back to Step 1 and compare the predicted $T_i^{pre}$ with the observed value $T^{obs}$. If the difference is larger than the tolerance, go to Step 2 in order to determine a new $B_j^{est}$. If the difference is smaller than the defined tolerance, provide the current $B_j^{est}$ as output.

Based on this algorithm, the derived travel times from pumping tests can be inverted using our self-developed code [41].

2.3. Wavelet De-Noising

A wavelet is a mathematical function used to separate a given function or continuous-time signal into different scale components [42,43]. Usually one can assign a frequency range to each scale component. Each scale component can then be studied with a resolution that matches its scale. A continuous wavelet transform (CWT) is the representation of a function by wavelets as Equation (7) [43]:

$$W_f\left(a,b\right)=\langle f,{\phi }_{a,b}\left(t\right)\rangle =\frac{1}{\sqrt{a}}\int f\left(t\right)\overline{\phi }\left(\frac{t-b}{a}\right)dt,$$

(7)

where $f\left(t\right)$ is the original signal; a is the scaling factor; b represents the time shift factor; and ${\phi }_{a,b}$ is a finite-length or fast-decaying oscillating waveform, which is known as “mother wavelet”.

Discrete wavelet transform (DWT) is an implementation of the wavelet transform using mutually orthogonal sets of wavelets defined by carefully chosen scaling and shift factors. It leads to an efficient iterative scheme for the transformation. Mallat [34] proposed a wavelet decomposition and reconstruction algorithm, which is known as a two-channel sub-band coder using quadrature mirror filters (QMFs) [44,45]. For an N-sample noised signal a0, the Mallat decomposition algorithm uses orthogonal filters of length L, and can be described as follows:

$$\begin{gathered} a_j\left(i\right)={\sum }_{m=0}^{L-1}{a}_{j-1}\left(2i-m\right)h\left(m\right), \\ {d}_j\left(i\right)={\sum }_{m=0}^{L-1}{a}_{j-1}\left(2i-m\right)g\left(m\right), \\ 0\le \mathrm{i}\le N_j-1, 1\le j\le \log N, \end{gathered}$$

(8)

where $a_j\left(i\right)$, $d_j\left(i\right)$ are the ith approximation and detail coefficient at level j, respectively. $h()$, $g()$ are low pass and high pass L-tap filters obtained from the chosen wavelet, respectively. $N_j=\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$2^j$}\right.$ is the number of wavelet coefficients in level j.

The Mallat reconstruction algorithm is [45]:

$$\begin{array}{c}{a}_j\left(i\right)={\displaystyle \sum _{n=i/2}^{\left(L-1+i\right)/2}}{a}_{j+1}\left(n\right)h\left(2n-i\right)+{\displaystyle \sum _{n=i/2}^{\left(L-1+i\right)/2}}{d}_{j+1}\left(n\right)g\left(2n-i\right)\\ 0\le \mathrm{j}\le \log N-1;0\le i\le N_j-1.\end{array}$$

(9)

Wavelet transforming and de-noising, following the Mallat algorithm, were applied in this study using the MATLAB^® wavelet analyzer app (Wavelet Toolbox Guide, Mathworks Inc.). The steps of the wavelet de-noising algorithm are given in Figure 2.

In order to remove the noise in a given target signal, Equation (8) is applied to decompose this signal into groups of coefficients at different frequency levels. Based on the different features of the target signal and noise in different frequency segments, certain threshold values are determined and applied to eliminate the noise. Thresholding can be applied by two different methods: hard and soft thresholding. In hard thresholding, elements are set to zero when their absolute values are lower than the threshold value and unchanged when higher than the threshold value. In soft thresholding, the coefficients that exceed the threshold value are softened and decreased by the threshold value. In this study we applied the soft thresholding method since soft thresholding provides smoother de-noising results, which meets the data requirement of hydraulic travel time inversion. After thresholding, the filtered coefficient sets are recomposed to obtain the de-noised signal based on Equation (9).

2.4. Diagnostic Plot Analysis

A diagnostic plot is a simultaneous plot of the pumping test induced groundwater level drawdown and the logarithmic derivative of the drawdown as a function of time in a log-log scale [39]. The mathematic expression of the logarithmic derivative of the drawdown as a function of time is shown as Equation (10):

$$\frac{\partial \mathrm{s}}{\partial \ln \mathrm{t}}=\mathrm{t}\frac{\partial \mathrm{s}}{\partial \mathrm{t}},$$

(10)

where s is the drawdown in the test well and t is the elapsed time since start of pumping.

The logarithmic derivative is highly sensitive to subtle variation in the shape of the drawdown curve and it allows for the detection of aquifer characteristics that are difficult to observe based on the drawdown curve alone. Therefore, the diagnostic plot (the joint use of the drawdown curve and its logarithmic derivative) is applied to characterize the aquifer and flow pattern geometry [46,47]. At the beginning of the plot, a unit slope indicates wellbore storage theoretically. A constant derivative at the immediate time represents a confined homogenous aquifer, while a U-shaped profile during intermediate times might represent an unconfined aquifer or a double porosity aquifer. A decreasing derivative in the late time plot implies a recharge boundary or a 3D flow geometry [46,47], while an increasing derivative indicates a no-flow boundary condition [39].

3. Field Work

3.1. Test Site Description

The test site is located in the north of Goettingen, Germany. It is at the eastern shoulder of the Leinetalgraben, where a number of folds, faults, and fractures are produced due to polyphase tectonic development under various tension forces and pronounced tectonics in the stories [48]. The test site comprises five 3” (ID 7.62 cm) PVC (Polyvinyl Chloride) wells that are arranged in a cross-shape configuration. The five wells are labeled N, S, W, O, and M for north, south, west, east, and middle, respectively, which represent the cardinal directions of the five wells. Each well is approximately 78 m deep and all five are constructed in the same way, with nine fully cased sections of 3 m length and nine screened sections of 5 m length in alternating order [49]. The screened sections are separated by clay seals within the well gravel pack to allow multilevel measurements. The well network is shown in Figure 3. The main rock types within the research area are limestone, siltstone, and claystone. Detailed well profiles for well O and well M are shown in Appendix A [50].

3.2. Multi-Level Pumping Tests Design

A series of multi-level pumping tests were conducted at the test site in April 2018. Figure 4 illustrates the experimental setup. The pumping tests followed a cross-well pattern, which means one well served as pumping (source) well and another well served as an observation (receiver) well. Groundwater was pumped out of well O using a submersible pump (Grundfos MP1, Bjerringbro, Denmark). This pump was connected to a frequency inverter in order to adjust the pumping rate. Pressure transducers were installed in both the pumping well (well O) and the observation well (well M) to monitor the drawdown response during the pumping tests. The drawdown data were recorded by a datalogger Campbell Scientific Micrologger CR3000 (Campbell Scientific, Logan, USA), which was connected to the pressure transducers. Double packer systems were used in both wells to allow depth specific pumping and drawdown measurements. Pressurized air was pumped into the packers to locate the pumping and measuring equipment at a specific depth and to prevent vertical flow within the well.

The double packer systems used in the pumping and observation wells are illustrated in Figure 5.

Following a tomographic configuration, the multi-level pumping tests were conducted in a sequential manner by positioning the pump packer system at a certain well screen and observing subsequently at multiple levels within the observation well. The pump did not work at shallow depth (i.e., at the first screen of the pumping well O) because of the relatively low hydraulic conductivity. Therefore, the pumping tests were conducted starting from the second screen to the fifth well screen. The second to the fifth screens of pumping well O were acting as source points, while the second to the fifth screens of observation well M were acting as receiving points in the tests (compare Figure 1). Given a 4 by 4 distribution in the sequential cross-hole configuration, 16 source-receiver pairs of data, with which 16 travel time data can be calculated, were obtained during the multi-level pumping tests. Figure 6 depicts the multi-level cross-hole array between the well O and the well M. The straight lines in the figure are not the real hydraulic signal travel paths. They demonstrate the spatial configuration of the source and receiver locations.

3.3. D Fracture Geometry Model

A 3D aquifer fracture geometry model of this test site was set up by Werner [51]. This model is based on geophysical data, distribution of calcite veins, paleostress field analysis, and microthermometry. Figure 7 shows an example of the rock samples which were obtained at the test site during the drilling process. With these data the possible positions and distributions of the fractures in the aquifer can be digitalized and visualized. Based on the result of Werner’s study, we set up a 3D fracture geometry model using the COMSOL^® Multiphysics software (COMSOL, Inc., Burlington, USA). An illustration of the 3D fracture geometry model is shown in Appendix B.

The inversion result obtained within our study will be compared with this 3D fracture geometry model.

4. Numerical Study

4.1. Numerical Modelling

A 3D numerical groundwater model, based on field data from the test site, was built to simulate the pumping tests performed in well O using the COMSOL^® Multiphysics software with the physics module “Darcy’s law”. The size of the investigated subsurface zone was set to 1.9 m × 1 m × 6 m, where 6 m is the length of one screen in the well O and 1.9 m is the distance between well O and well M. A mass flux point was set in the pumping well to simulate the pumping process and an observation point was set at the horizontally corresponding position of the observation well (1.9 m away from the pumping well). The mass flux was set to 10.2 L/min, according to the average pumping rate in the field experiments. Infinite element domains were established in all directions outside the investigated zone in order to eliminate boundary effects. An illustration of the numerical model setup is shown in Figure 8. The one-screen investigation zone model is a simplification of a multi-level pumping test at the test site, of which the depth orientated pumping and the water head measurement are at the same depth, with the assumption that the double-packer systems can totally prevent vertical flow within both pumping and observation wells. Thus, the remaining part of the well does not contribute to the pumping test. The aim of this modelling is to find the optimal de-noising parameters. The 3D-effects on the head response due to the other screens are not considered in this simulation.

The initial setting of the hydraulic parameters in the model was based on the former investigation work at the test site [52,53]. In order to ensure that the numerical model provides a similar hydraulic behavior compared to the real test site, the hydraulic parameters of the model were calibrated to match the drawdown and the calculated diagnostic plot to the observed data and diagnostic plot from the field experiment (Figure 9). The final setting of the model parameters is shown in

Table 1. Parameters of the numerical model.

Parameter		Value
Well skin	Porosity	0.4
	Hydraulic conductivity (m/s)	$1.7\times {10}^{-5}$
	Specific storage (1/m)	$5\times {10}^{-4}$
Matrix	Porosity	0.2
	Hydraulic conductivity (m/s)	$2\times {10}^{-8}$
	Specific storage (1/m)	$3\times {10}^{-3}$
Fracture	Porosity	0.4
	Hydraulic conductivity (m/s)	$8\times {10}^{-5}$
	Specific storage (1/m)	$1\times {10}^{-5}$
	Fracture thickness (m)	0.01

. The diagnostic plots of one set of real experimental drawdown data from pumping well O (Figure 9a) and the simulated data from the numerical model (Figure 9b) are shown below. The black curve is the drawdown curve while the blue curve shows its derivative with respect to time.

As shown in Figure 9, the logarithmic derivatives in both plots increase rapidly with a slope of approximately 1 from 1 s to 8 s, which indicates the presence of borehole storage effects. The U-shaped profiles of both derivative curves after 10 s represents double-porosity properties, which is consistent with the porous-fractured aquifer character of the test site. In general, the diagnostic plot of the numerically simulated drawdown data shows the same pattern as that from the field data, which indicates that the model aquifer has similar properties as the real one. Thus, the simulated data from the numerical model can be used to find a suitable method to de-noise the field experimental data.

In the model, an observation point corresponding to the position of well M was also set, which is 1.9 m away from the pumping well O. The simulated drawdown response at the observation point is shown in Figure 10a. Figure 10b shows the first derivative of the drawdown with respect to time. As shown in Figure 10b, the time t_peak at the derivative maximum, a measure for the hydraulic signal travel time, is 14.6 s.

4.2. Comparison Between Wavelet and Other Filters

Gaussian white noise was added to the numerically simulated drawdown curves in order to compare the results of the wavelet de-noising method with six filters provided with the MATLAB^® (The MathWorks, Inc., Natick, USA) “smooth” command (the names of those filters are moving, lowess, loess, sgolay, rlowess and rloess).

Wavelet de-noising was performed with the integrated “Wavelet Analyzer” of MATLAB^®. Mother wavelet “Symlets” was selected and the mother wavelet level was set to eight. The de-noising level was set to five. Figure 11 shows the original signal with a 30 signal-to-noise ratio Gaussian white noise, the de-noising result of wavelet de-noising, and the de-noising results of the six different filters.

The travel time value t_peak of each de-noised signal was calculated. The calculation results are shown in

Table 2. Calculated travel time values t_peak.

Filter	tpeak (s)
Wavelet	15.2
Moving	15.8
Lowess	15.8
Loess	13.2
Sgolay	15.8
Rlowess	15.8
Rloess	28.6

.

As shown in the table, the t_peak value derived from the result of wavelet de-noising is the closest to the known “true” t_peak value, which is 14.6 s (Section 4.1). The difference among the results of the seven de-noising methods is not significant. However, since the t_peak value in real experiments can be very small, small differences can cause large deviations of final results. Therefore, although the lead is not prominent, we conclude that the wavelet de-noising method has the best de-noising performance and that it is suitable to be utilized to de-noise the experimental data obtained from the pumping tests.

4.3. Wavelet De-Noising Method

So far, it is shown that the wavelet de-noising method is suitable to de-noise the experimental drawdown data obtained from the pumping tests. However, which mother wavelet and which de-noising level should be selected to remove the noise in the experimental drawdown signal and to get the most accurate travel time value is not yet known. The simulated drawdown response and the already known time t_peak from the numerical model were used to find the best wavelet de-noising method to de-noise the experimental data. Three variables were investigated: wavelet type, mother wavelet level, and de-noising level.

Since the experimental conditions may change from day to day, the noise levels in the drawdown data can be different. Therefore, Gaussian white noise with different signal-to-noise ratios (10, 20, 30, and 40) was added to the simulated drawdown response from the numerical model.

Different mother wavelet types were tested in this study and four different types of mother wavelets were found effective for our work. The names of these four mother wavelets are “Symlets”, “Coiflets”, “Biorthogonal”, and “Dmey” (Wavelet Toolbox Guide, The MathWorks Inc., Natick, USA).

The level of the mother wavelet is a parameter of the de-noising method. In order to find the influence of the mother wavelet level on the de-noising result, the simulated drawdown curve with a signal-to-noise ratio of 20 noise was selected as the de-noising object, “Symlets” was selected as the mother wavelet, and the de-noising level was set to five, the mother wavelet level was varied from four to eight to de-noise the signal (level eight is the optional maximum level in MATLAB^® wavelet app), and the de-noising results are shown in Figure 12.

As shown in Figure 12, there is not a particularly big difference in the de-noising results with different mother wavelet levels. Therefore, the mother wavelet level is not a sensitive parameter in this study. It also holds if the de-noising level is changed. However, if we zoom into the figure and focus on the beginning part of the curve (Figure 13), which is the most important part for travel time derivation, it is found in the figure that the higher the mother wavelet level, the smoother the signal curve after de-noising. Since the drawdown curve should be a smooth curve, and the first derivative is required for seeking the travel time, the level of the mother wavelet level was selected to the highest level in the following study.

The type of the mother wavelet and the de-noising level are the two other parameters in the de-noising method. Therefore, the drawdown data with different noise levels (signal-noise-ratio (SNR) = 40, 30, 20, 10) were de-noised with four different mother wavelets, and with different de-noising levels (level = 5, 6, and 7) respectively. The de-noising results were then used to calculate the travel time t_peak using MATLAB^®. The de-noising results with lower de-noising levels (i.e., lower than level 4) were not smooth enough and can hardly be used for further calculation of the t_peak.

The results are summarized in

Table 3. Travel times t_peak for different mother wavelet types and de-noising levels.

SNR = 40
De-Noising Level	Mother Wavelet
	Symlets	Coiflets	Biorthogonal	Dmey
Level 5	12.4 s	19 s	23.4 s	15 s
Level 6	22.4 s	15.2 s	22.8 s	17.2 s
Level 7	22.4 s	23.6 s	22.4 s	27.6 s
SNR = 30
De-Noising Level	Mother Wavelet
	Symlets	Coiflets	Biorthogonal	Dmey
Level 5	12.8 s	13.4 s	13 s	13.8 s
Level 6	23.8 s	15.4 s	23.8 s	18 s
Level 7	22.8 s	24 s	23.2 s	27.6 s
SNR = 20
De-Noising Level	Mother Wavelet
	Symlets	Coiflets	Biorthogonal	Dmey
Level 5	/	/	/	/
Level 6	25 s	13.2 s	23.4 s	17.6 s
Level 7	21.6 s	22.2 s	22 s	26.8 s
SNR = 10
De-Noising Level	Mother Wavelet
	Symlets	Coiflets	Biorthogonal	Dmey
Level 5	/	/	/	/
Level 6	/	/	/	/
Level 7	21.2 s	20.6 s	21.8 s	24.4 s

. Based on the numerical model simulation, the exact travel time t_peak should be 14.6 s (Section 4.1).

The slash in the table means that no unique t_peak can be found from the first derivative curve, because of the high level of noise and poor de-noising result.

After analyzing the results in Table 3, several conclusions can be drawn:

In general, the results using mother wavelets “Coiflets” and “Dmey” are relatively better than the ones using “Symlets” and “Biorthogonal”.

Mother wavelet “Dmey” is suitable for de-noising less noisy signals with low de-noising levels (level five particularly).

Mother wavelet “Coiflets” behaves the best when the noise is relatively large and the de-noising level is six.

When the de-noising level is seven, the relative error of the derived travel time t_peak is relatively large, i.e., the derived travel time value is not accurate.

When the signal-to-noise ratio is smaller than 20, wavelet de-noising cannot effectively remove noise signals. Therefore, the travel time value calculated from the de-noising result is not reliable.

Although the noise is added artificially and the noise is merely Gaussian white noise, the wavelet de-nosing method used in this study is not capable to completely remove the noise and to get the exact ‘true’ travel time value.

In the actual field drawdown data, the noise is generally large and the signal-to-noise ratio is usually between 20 and 30. It can be found from Table 3 that whether the signal-to-noise ratio is 20 or 30, the travel times derived using “Coiflets” mother wavelet and de-noising level six are closest to the ‘true’ travel time value. Therefore, mother wavelet “Coiflets”, with mother wavelet level five and de-noising level 6 are used in this study to process the field data.

5. Results

5.1. Travel Time Results from Experimental Data

Using the wavelet de-noising parameters defined above, the field data from multi-level pumping tests were processed to remove the noise signal, and the de-noising results were used to obtain the travel time values t_peak. Figure 14 shows an example of data processing of a field data set. The t_peak results of all pumping tests are summarized in Figure 15. The figure also shows the calculated distance between the pumping source and the observation point based on the well construction profiles and well arrangement (in bracket). The travel time results for O2M5 (the second screen of well O as a source and the fifth screen of well M as a receiver), O3M5 (the third screen of well O as a source and the fifth screen of well M as a receiver), and O4M5 (the fourth screen of well O as a source and the fifth screen of well M as a receiver) are missing, because the noise of these field data sets is too high to obtain reliable travel time values.

As shown in Figure 15, since the distance between the signal source and receiver location is relatively small, the signal travel time for horizontal directions (e.g., O2M2) is relatively short. The signal travel time for experiments with a high slope between the signal source and the receiver locations (e.g., O2M4) is relatively large because of the longer travel distance. However, the difference in travel time can be observed in the case of a similar distance between the source and receiver. For instance, the travel time from the fourth screen of well O to the second screen of M (6.4 s) is significantly smaller than the travel time from the second screen of well O to the fourth screen of M (25.6 s), while the travel distances within the two experiments are basically identical. The explanation is that the hydraulic connection between the fourth screen of well O and the second screen of M is better compared with the other setup, due to a possibly existing fracture.

5.2. Hydraulic Travel Time Inversion Result and 3D Fracture Geometry Model

The diffusivity distribution based on inversion using the SIRT algorithm is shown in Figure 16. The main features of the inversion result are the three-horizontal high-diffusivity channels between well O and well M. The diffusivity of the cross section is relatively high near each well screen.

Figure 17a shows the result of a geological 3D fracture geometry model of the top 45 m based on the data introduced in Section 3.3 and Figure 17b shows the 2D travel-time based inversion result for comparison purposes.

The grey lines in Figure 17b are possible locations of fractures based on the 3D geometry model. In general, the zones where a fracture may exist according to the 3D geometry model have higher diffusivity values in the inversion result. Although there are slight deviations, the generally consistent fracture locations and high diffusivity zones indicate a reliable diffusivity distribution within the study area obtained from the travel-time based hydraulic tomography.

6. Conclusions

In this study, travel-time based hydraulic tomography was used to identify the two-dimensional hydraulic diffusivity distribution between a pumping and a drawdown observation well at our test site. The agreement between the obtained hydraulic diffusivity distribution and the geological 3D geometry model proves the reliability of the hydraulic tomography result. The employed travel-time based inversion requires accurate hydraulic signal travel time values, which determine the quality of the inversion result.

The major obstacle in getting accurate travel times is the noise in the experimental drawdown curves. In order to find the best way to de-noise the drawdown curves, a numerical 3D model based on our test site data was established in this study. Using diagnostic plots, it was demonstrated that the model provides a hydrogeological behavior according to the test site situation. By using wavelet and other filters to de-noise numerically simulated pumping test drawdown curves, with added Gaussian white noise, the wavelet de-noising method, with mother wavelet “Coiflets”, mother wavelet level five, and de-noising level six, proved to have the best performance in processing the noisy data.

Although the best de-noising method was found in this study, it cannot totally remove the noise, which might still influence the accuracy of the result. In addition, this study assumes that the noise in the experimental data was white noise. If the experimental data contain other types of noise, the de-noising method introduced in this study should be further tested and if needed improved.