Freely available online through the SEG open-access option.

ABSTRACT

Applications of time-lapse inversion of electrical resistivity tomography allow monitoring variations in the subsurface that play a key role in a variety of contexts. The inversion of time-lapse data provides successive images of the subsurface properties showing the medium evolution. Image quality is highly dependent on the data weighting determined from the data error estimates. However, the quantification of errors in the inversion of time-lapse data has not yet been addressed. We have developed a methodology for the quantification of time-lapse data error based on the analysis of the discrepancy between normal and reciprocal readings acquired at different times. We applied the method to field monitoring data sets collected during the injection of heated water in a shallow aquifer. We tested different error models to indicate that the use of an appropriate time-lapse data error estimate yielded significant improvements in terms of imaging. An adapted inversion weighting for time-lapse data implies that the procedure does not allow an over-fitting of the data, so the presence of artifacts in the resulting images is greatly reduced. Our results determined that a proper estimate of time-lapse data error is mandatory for weighting optimally the inversion to obtain images that best reflect the evolution of medium properties over time.

INTRODUCTION

Time-lapse electrical resistivity tomography (ERT) has become a widely used method to characterize hydrogeologic processes. The method has been used to study water and snowmelt infiltration in soil and in the vadose zone (e.g., Daily et al., 1992; LaBrecque et al., 2002; French and Binley, 2004; Oberdörster et al., 2010; Slater et al., 2010; Travelletti et al., 2012), groundwater flow (e.g., Coscia et al., 2011; Doetsch et al., 2012), solute transport (e.g., Kemna et al., 2002; Koestel et al., 2008; Ogilvy et al., 2009; Robert et al., 2012; Rucker 2014), changes in the thermal state of permafrost-affected rocks (e.g., Krautblatter et al., 2010), heat conduction and transport (e.g., Hermans et al., 2015), and many others.

It is well-known that ERT suffers from loss of resolution with increasing distance to the electrodes (e.g., Nguyen et al., 2009; Perri et al., 2012), which can be analyzed using sensitivity and uncertainty analyses (see, e.g., Hermans et al., 2012b; Robert et al., 2012). Furthermore, variations in the data error during monitoring (e.g., variations in the galvanic contact of the electrodes) may also lead to electrical images of variable quality in time, potentially misleading their quantitative interpretation (Deceuster et al., 2013). Therefore, several approaches have been suggested to improve the reliability of the time-lapse images by including a temporal regularization (e.g., Kim et al., 2009; Hayley et al., 2011; Karaoulis et al., 2011; Loke et al., 2014). Such regularization constrains the inversion to a reduced range of acceptable models, increasing the correlation of inverted parameter values with hydrogeologic parameters of interest (Loke et al., 2014). However, to date, no study has focused on the impact of the time-lapse data error description in the inversion, although it has been recognized for being of critical importance (e.g., Miller et al., 2008; Flores-Orozco et al., 2012). The objective of this paper is to develop a methodology for estimating the time-lapse data error applied for the data weighting in the inversion.

The comparison of normal and reciprocal measurements is commonly used to quantify the data error in static imagery (e.g., LaBrecque et al., 1996; Slater et al., 2000; Koestel et al., 2008; Oberdörster et al., 2010; Flores-Orozco et al., 2012). The reciprocal measurements correspond to a swap of the electrodes used for current injection and voltage measurements during the “normal” acquisition. Here, we assume that the reciprocal theorem applies on the resistance differences and we analyze normal-reciprocal data-difference discrepancies to estimate appropriately the weighting of time-lapse ERT data in the inversion. We investigate the impact of the data error weighting on the reconstructed images. Our results demonstrate that a specific data error estimation method is required for time-lapse inversions to obtain reliable images in ERT monitoring.

ERROR IN TIME-LAPSE ERT

The error in ERT measurements can be considered as being composed of a random component and a systematic component correlated over time. At a time t0, data d0 can be estimated from the model parameters m0 by the computation of the nonlinear function f(m0):  
d0=f(m0)+εS+ε0.
(1)
In the inversion, log10-transformed resistances are used as data d and log10-transformed resistivity (of lumped finite-element cells) as parameters m to account for the typically large range of resistivity values for earth materials. The functions εS and ε0 stand for the systematic and random error at t0, respectively (LaBrecque and Yang, 2001). Random error cannot be predicted because it arises primarily from fluctuations in the contact between the electrodes with the ground/air, and in the injected current and its pathways (e.g., Binley et al., 1995; LaBrecque et al., 1996; Slater et al., 2000). Systematic error might be associated with problems during data acquisition such as poor galvanic contact, malfunction of the measuring device, misplaced, or misconnected electrodes; but it also might be due to numerical errors. In the case of monitoring, further sources of error are related to changes in the position of the electrodes and their contact with the ground for measurements collected at different times due to the medium saturation variations, for instance, also associated with variations in the injected current and signal strength (e.g., Supper et al., 2014).
For time-lapse imaging, LaBrecque and Yang (2001) demonstrate that the data-differences inversion method allows reducing the effect of systematic error present in the data, which cancels out by subtraction. ERT measurements typically present a better repeatability compared with the data value accuracy, expressly when electrodes are set up at the beginning of the experiment and remain in place during monitoring (LaBrecque and Yang, 2001). In the case in which the correlated error dominates, it is particularly worth performing data-differences inversions seeking changes in the medium from the initial state. At time ti, the data differences from the initial state are written  
did0=f(mi)f(m0)+εiε0,
(2)
where di and mi represent, respectively, the measured data and the model parameters at time ti. Because the uncorrelated error εi at time ti is independent of ε0, the error variance of the data differences corresponds to the quadratic addition of the respective variance of ε0 and εi. Therefore, data-difference inversions are of no interest for systems presenting uncorrelated source of error with a higher amplitude than the systematic error.

The difference inversion algorithm works in two steps: First, a background model m0 is reconstructed from data acquired at the initial state t0. The error affecting the background image is due to the propagation of the random and the systematic error component. Usually the data error is deduced by comparing normal and reciprocal measurements (LaBrecque et al., 1996; Slater et al., 2000; Koestel et al., 2008). During a second step, the parameter changes mim0 are sought by inverting the data differences did0. The difficulty of this approach is to estimate appropriately the uncorrelated error level because the fractions of the random and systematic error are rarely known (LaBrecque et al., 1996). Furthermore, classical data error estimation procedures using reciprocal measurements are not adapted to the time-lapse measurement scheme. For instance, data with large geometric factors are characterized by low signal-to-noise ratio and often associated with larger weights because the classical error model applied for background inversion presents an increase of the error level with increasing resistances (Slater et al., 2000). Hence, the difference inversion with classic error estimation procedure typically leads a reduced influence of such measurements, or a severe smoothing in the concerned sensitive regions. To improve the difference inversion scheme, an adequate fitting to the data differences and their associated error needs to be performed for weighting adequately the data in the inversion. To our knowledge, there exists no study dealing with the quantification of error associated with differences of time-lapse data. Currently, the standard procedure corresponds to the static error assessment described below (Slater et al., 2000; Koestel et al., 2008). However, the static error model does not take into account that data differences are actually fitted in the inversion, as performed here, to remove the correlated error affecting the data. Hermans et al. (2016) stress the importance of defining a specific measurement error for time-lapse acquisitions from time-lapse resistance changes observations. However, their study is based on synthetic data only and they use a prediction-focused approach, so no inversions are carried out.

Here, we adopt the common understanding that the normal-reciprocal discrepancy can be considered as a practical measure of data error (e.g., LaBrecque et al., 1996; Slater et al., 2000; Flores-Orozco et al., 2012). A linear error model for the resistance has been proposed to describe the static data error εB=εS+ε0 used for the background image reconstruction (LaBrecque et al., 1996). That resistance error model is given by (Slater et al., 2000)  
εB=a+bR.
(3)
Following the description from Slater et al. (2000), in our study, we define the parameters a and b, such that the error model defined in equation 3 encompasses all data points in a logΔR versus logR¯ plot (Figure 1). The value ΔR corresponds to the difference between resistances measured in a normal and a reciprocal configuration, whereas R represents the mean value of normal and reciprocal resistances. Considering that the error model determined in such a way is mainly defined by the larger normal-reciprocal misfits, its choice implies a conservative fit of the data to prevent the creation of artifacts.

MATERIALS AND METHODS

Time-lapse error modeling

Here, we propose to extend the analysis of misfit between normal and reciprocal for data differences, i.e., the data given by the difference between normal background measurements and those collected at ti and their corresponding reciprocals. As for the static case, we examined how the normal and reciprocal measurements variations with time could be used for quantifying the error. In equation 2, di and d0 are expressed in log10 resistances, so the time-lapse error should correspond as well to a difference of log10 resistances. Because we perform a data-difference inversion, we need to characterize how the normal and reciprocal measurements evolve with time from the initial acquisition. The comparison of the time-lapse measurement variations |logRN,ilogRR,i| versus the initial state variations |logRN,0logRR,0| is not adapted for estimating the time-lapse error as illustrated by the wide dispersal of the data on the scatterplot in Figure 2a. Instead, the scatterplot showing the reciprocal resistances differences ΔlogRR, i.e., the quantity that is inverted in the data difference inversion, as a function of the normal resistances differences ΔlogRN shows a better 1:1 consistency (Figure 2b). The value ΔlogRN corresponds to the normal (N) difference between t0 and ti, and ΔlogRR represents the reciprocal (R) difference, such that  
ΔlogRN=logRN,ilogRN,0,ΔlogRR=logRR,ilogRR,0.
(4)
Contrasted behaviors are observed on the ΔlogRR versus ΔlogRN scatterplot related to the different phases of the experimental setup as described below. Nevertheless, the observed correlation between ΔlogRN and ΔlogRR shows that the error quantification for time-lapse inversion can be deduced from their comparison. Indeed, the plot of the time-lapse normal-reciprocal misfit, i.e., the difference between ΔlogRN and ΔlogRR, reveals a clear increase in the time-lapse data error with decreasing resistance values (Figure 3). Thus, we suggest that the data-differences error εTL,i=ε0+εi, i.e., the normal-reciprocal discrepancy, can be modeled as  
|ΔlogRNΔlogRR|=aR+b.
(5)
We propose here to fit the parameters a and b in equation 5 in three different strategies:
  1. 1)

    As an envelope fit, such that the error model encompasses all data points in a |ΔlogRNΔlogRR| versus logR¯i plot, implying a conservative fit of the data similar to the above-mentioned approach of Slater et al. (2000). This is achieved by computing the mean values and standard deviation that |ΔlogRNΔlogRR| takes in bins (one bin per decade of resistance values) and fitting equation 5 in the least-squares sense to the means plus two standard deviations (encompassing 97% of the data).

  2. 2)

    A least-squares fit to the whole data set.

  3. 3)

    A resistance independent error model (a=0 in equation 5).

The results of the inversions using the three types of time-lapse error models are compared, and the reliability of the different obtained images is discussed later.

Inversion algorithm

The inversion is performed with CRTomo (Kemna, 2000) that uses a Gauss-Newton scheme to iteratively minimize the objective function Ψ(m):  
Ψ(m)=Ψd(d,m)+λΨm(m),
(6)
where Ψd(d,m) is a measure of the data misfit, the model functional Ψm(m) defines desired model constraints, and λ is the so-called regularization parameter. We refer to Kemna (2000) for a complete description of the inversion algorithm.
The data misfit Ψd(d,m) expresses differently for the reconstruction of the background or the time-lapse images. For the static case, we assume that each resistance measurement dj is contaminated by an error εB,j, so that the measure of the error-weighted data misfit for N data points is defined as  
Ψd(d,m)=j=1N|djfj(m)|2|εB,j|2.
(7)
For the time-lapse inversion, the measure of the data misfit Ψd(d,m) is defined using a double difference as stated in equation 2. The first difference compares the background data d0 and the monitoring data di, whereas the second one evaluates the discrepancy inferred from the forward model estimations computed with the background state model f(m0) and the monitored model f(mi):  
Ψd(d,m)=j=1N|(di,jd0,j)(fj(mi)fj(m0))|2|εTL,j|2,
(8)
with εTL,j, the data-differences error estimated as described in equation 5. Such a definition of the data misfit Ψd(d,m) has been demonstrated to improve imaging results for the inversion of time-lapse ERT data sets (e.g., LaBrecque and Yang, 2001; Kemna et al., 2002).
The model functional Ψm(m) allows the spatiotemporal regularization of the inversion  
Ψm=Wm(mm0)2.
(9)
To focus solely on the data weighting impact, the regularization term contains no explicit temporal constraint (i.e., no additional parameter as in Kim et al., 2009; Karaoulis et al., 2011; Loke et al., 2014). For the background reconstruction m0=0, whereas for the time-lapse inversion, m0 corresponds to the background model (e.g., Twomey, 1963) inducing therefore an implicit temporal constraint.
During the inversion, a univariate line search is performed at each Gauss-Newton iteration step to find the optimum value of the regularization parameter λ that balances the respective weights of the data reduction Ψd(d,m) and the model smoothing Ψm(m). The line search seeks the λ value that minimizes the residuals σ that yields the desired target misfit: σ=|1ξ|. We used here a ξ value of 102. For the background inversion residuals σ0, we write  
σ0=1Nj=1N|djfj(m)|2|εB,j|2,
(10)
whereas the time-lapse inversions residuals σTL read as  
σTL=1Nj=1N|(di,jd0,j)(fj(mi)fj(m0))|2|εTL,j|2.
(11)
Such a definition of the inversion stopping criteria avoids over-fitting the data and thus the presence of artifacts in the resulting images. The explicit formulation of a final misfit threshold that constrains the dynamic adjustment of λ is explicitly defined in Occam’s inversion method (Constable et al., 1987; Aster et al., 2013). Therefore, the data error estimation is of particular importance because it affects not only the data misfit weighting but also the satisfactory level of misfit at which the inversion is stopped.

Experimental set up

ERT measurements were collected during a heated water injection in sands at shallow depth (5 m). Detailed descriptions of the study area are given in Vandenbohede et al. (2011) and Hermans et al. (2012b). The description of the experiment itself can be found in Hermans et al. (2012a). The site is located on the campus De Sterre of Ghent University, Belgium (Figure 4). The schematic cross section of the site is presented in Figure 4 and shows sands from 0 to 4.4  m and a sandy clay layer below. The water table is nearly flat at the site at 2  m, with a hydraulic gradient of 0.005 toward the south–southwest, as derived from measurements in three wells (W01, W02, and W03). The temperature measured in the aquifer during the ERT monitoring was 10.5°C (Hermans et al., 2012b).

Injection of tap water and heated water was performed on the site, and surface ERT measurements were collected before, and during the injection to monitor the flow and transport of heat (associated with the heated water) in the subsurface. Tap water was injected initially over 15 days to avoid mixing cold groundwater and heated tap water. Then, heated tap water (49°C) was injected for 144 h (six days, starting on 18 January 2011) at a rate of 80 l/h and at a depth of 4 m. The ERT monitoring line was centered on the injection well with 48 electrodes spaced by 60 cm (28.2 m long profile). The cables and the electrodes were left in place during the heating phase. A dipole-dipole configuration was used with a dipole length of a=3 and a separation between current and potential dipoles n13, leading to 579 measurement points. Data were collected as normal and reciprocal pairs by setting accordingly measuring protocols with an ABEM SAS1000 system using a constant injection current of 200 mA. The duration of a whole data set acquisition was of 3 h including normal and reciprocal measurements.

RESULTS

Data analysis

The pseudosection for the background is shown in Figure 5a. The pseudosection shows higher apparent resistivity values closer to the surface and a monotonic decrease at larger depths. For data collected at different times, we present the pseudosections of percentage changes of the apparent resistivity values. The percentage change was computed for each quadrupole as the change in apparent resistivity at each time-lapse ti with respect to the background readings. Pseudosections of percentage changes show large variations and no obvious systematic error. Note that two effects play a role here. The first one is the injection of the heated water, which produced localized changes around the borehole marked by a decrease of the apparent resistivity. The second one is related to a saturation decrease in the vadose zone. Before the background data acquisition, heavy rains occurred and water infiltrated the unsaturated zone. The shallow part of the medium was more saturated for that background acquisition than for the subsequent ones because the water could circulate easily through the shallow sandy medium toward the saturated zone. The process of saturation decrease induces an increase of the apparent resistivity and yields changes in the entire pseudosection given the integrative nature of the measurements.

Considering the expected low velocity of the diffusive process (Vandenbohede et al., 2011), we could collect reciprocal measurements for each time frame of the cooling phase. For the background error model, the estimated parameters’ values of equation 3 are a=1.8  mΩ and b=0.7 (Figure 1b). The amplitude of the constant parameter a=1.8  mΩ of the linear error model allows wrapping the data error estimated for low resistance values away from the 1:1 line on the scatterplot (Figure 1). In case the acquisition of the full reciprocal is not possible due to rapid changes in the subsurface, we suggest to use a subset of reciprocals representative of the whole data sets (e.g., Hermans et al., 2012b).

For the time-lapse noise characterization, the ΔlogRR versus ΔlogRN scatterplots present different behaviours depending on the acquisition time (Figure 2b). Data corresponding to the first acquisition during the heat-plume injection (t=24  h) present low variations from the background acquisition. The time-lapse data are closer to the initial state because less volume has been injected compared with the following acquisitions. Therefore, the difference values from the t0 acquisition are lower and present a dispersed distribution. The data variations from the initial state corresponding to the t=72  h acquisition present a more focused distribution, likely due to the development of the heat plume in the medium. For the acquisition at t=144  h, data differences present a higher range of variability and an alignment with the 1:1 line. For that time, data are affected by superficial processes and by the influence of the heat-plume development. The trend of the time-lapse reciprocal misfit has been fitted using equation 5 with different models shown in Figure 3. Thus, an envelope curve for which the parameters’ values are estimated to be a=0.001 and b=0.09 allows wrapping the whole data sets. For the least-squares fit, a=7×106 and b=0.1. In addition, a constant error model independent of the measured resistance was also selected with b=3. The least-squares and the envelope models present a similar value for the parameter b that corresponds to the minimal error for high resistance values. By the contrary, the parameter a presents consequent variations larger than two orders of magnitude when comparing the least-squares and the envelope models.

Imaging

The electrical resistivity image for data collected before the injection of heated water is presented in Figure 6. To better solve for the horizontal layering observed in the sedimentological/hydrologic units (Hermans et al., 2012b), the horizontal/vertical anisotropy ratio in the smoothness constraint inversion was set to 1000∶1 at 2 and 4.5  m depth, considering horizontal limits observed in borehole sediments. It means that the smoothing effect of the inversions is 1000 lower vertically than horizontally, permitting the reduction of the vertical smoothing effects at the limit between layers (Caterina et al., 2014). Imaging results, computed with such constraint, show agreement with the sedimentological/hydrologic layering because it can be observed in Figure 4.

The first approach used for the inversion of time-lapse data sets was to define the error as a constant error b=3, independent of the measured resistance values (Figure 7). The images show a conductive plume in the center, over-smoothed and hardly visible. We also see an increase of resistivity in the upper part of the medium that corresponds to the saturation decrease effect. We note that the regions of resistivity increase are blurred and seemed to penetrate deeper than the water table located at 2  m, where no saturation variations are expected.

The inversions performed using a resistance-dependent error model with a least-squares fit as defined in equation 5 (Figure 3) are shown in Figure 8, whereas the ones performed with an envelope fit are shown in Figure 9. The images shown in Figure 8 reveal several artifacts, which impede the interpretation of the results. Clearly, the inversions tend to over-fit the data. Note that for these inversions, the target data misfit is never reached. The use of the bins median (Koestel et al., 2008) instead of a least-squares estimate provides similar results in terms of imaging (not shown here).

The inversions shown in Figure 9 reveal the heat injection plume as well as the changes in the unsaturated layer. After 72 h, the minimum anomaly corresponding to the heat plume ranges from 10% to 20% with a continuous increase of the spatial extension of the heat plume at a depth of approximately 5  m. The heat plume reaches its maximum extension after 144 h, associated with a resistivity decrease of values up to 35%. In addition, regions affected by an increase of resistivity related to the decrease of the saturation in the vadose zone are limited to the shallow region, above the water table. The envelope error model thus produces images that also better reflect shallow phenomena. Hence, the results shown in Figure 9 help to demonstrate the relevance of an adequate error parameterization for the quantitative imaging of time-lapse ERT data sets.

The inversion results of time-lapse ERT are influenced by the weighting of the data misfit. A characterization of the error model performed from the analysis of the normal and reciprocal data variations adapted to time-lapse experiments is important for the objectiveness of the interpretation. Moreover, the estimation of the error perturbing the data in time lapse from the reciprocal data highlights the significance of performing such measurements despite the increasing time of acquisition. Further regularization schemes, such as those presented in Kim et al. (2009) or Karaoulis et al. (2011), may help to refine the interpretation. However, the choice of regularization method requires a priori information on the spatiotemporal evolution of the medium properties that are not always available. Hence, one challenge in the 4D methodology is that independent inversions are first run to estimate the value of the constraints that are subsequently used. This step is decisive to distinguish if the resistivity variations present a sharp or a smooth pattern because they require different regularization constraints (Loke et al., 2014). Therefore, the error level estimates at this very first step have a critical influence on the 4D constraint, whatever the inversion algorithm.

CONCLUSION

We performed a detailed analysis of time-lapse ERT measurements using reciprocal readings. In our study, we found a dependence of the time-lapse error on the magnitude of the measured resistance (or the measured voltage at a constant current). Based on our analysis, we proposed a methodology for the error model estimation to describe the data error present in the data differences between the baseline and time-lapse measurements. We performed inversions of time-lapse data using different error model parameterizations, using an envelope fit, a least-squares fit, and a constant error model. We found that significant improvements are obtained by inverting data differences based on a resistance-dependent error model using the envelope approach. The reduction of acquisition time is decisive when monitoring rapid changes of the medium properties, and one may consider using a limited number of reciprocal acquisitions.

We therefore conclude that the analysis of the normal-reciprocal misfit for time-lapse differences is a suitable tool for the evaluation of time-lapse data quality and to quantify data error. The characteristics of the measuring protocol may have a strong influence in the error parameterization, but the suggested error model is dependent on the values measured, like previous models for static measurements. Moreover, the error model estimate could be applied for measurements with high dynamics in the measured resistances or for configurations with a low dynamic in the measured values (where the error model tends to a constant value). Our results suggest that further improvements in the regularization approach, i.e., inclusion of a prior model, spatiotemporal constraints, etc., may also provide better results when the data error is adequately quantified. Indeed, the data need to be weighted, no matter the inversion procedure. Therefore, the findings in our study are relevant for all applications of ERT monitoring, independent of the inversion algorithm.

ACKNOWLEDGMENTS

The work by N. Lesparre was supported by the project SUITE4D from a BEcome a WAlloon REsearcher fellowship fund cofinanced by the Department of Research Programs of the Federation Wallonia — Brussels and the COFUND program of the European Union.