Time-lapse inversion is frequently used to map changes in electrical conductivity through time, as might be expected from processes such as infiltration of water into dry soils or tracer transport (Figure 14). Here, we demonstrate the steps taken for data collection in the field and the impacts of inversion decisions using a case study at the former Brandywine Defense Reutilization and Marketing Office (DRMO), which is an inactive 3-hectare facility administratively controlled by Andrews Air Force Base (AAFB) and located approximately 13 kilometers south-southeast of AAFB in Brandywine, Maryland, USA. The Brandywine DRMO yard was used from 1943 to 1987 for temporary storage of scrap materials and hazardous waste generated from various Department of Defense facilities in the region and is currently classified as a Superfund site. The primary groundwater contaminant at the former DRMO is trichloroethylene, which has spread beyond the AAFB property into adjacent commercial and residential properties (United States Air Force, 2006). We note that this case study shows the state of the art in terms of data collection, inversion, and analysis, especially for a 4-D system, rather than the state of the practice, which remains 2-D static inversions at the time of this writing.

Figure 14 – Schematic of a cross-well inversion from a system with active changes—in this case, the movement of an electrically conductive tracer. By looking at a series of inversions through time, we would see changes in conductivity at particular pixels, as demonstrated by the “pixel breakthrough curves” (e.g., Slater et al., 2000) on the right.

To remediate groundwater contaminants at the site, a biostimulation effort commenced in 2008. This effort involved injecting amendments into the subsurface to enhance microbial activity associated with contaminant biodegradation. An autonomous cross-well ER monitoring system was installed at the site to demonstrate and validate the utility of geophysical measurements for providing timely and actionable information on the spatiotemporal behavior of amendments injected into the subsurface.

ER Monitoring System

The Brandywine ER configuration consists of seven boreholes with 15 electrodes each for a total of 105 electrodes (see Figure 15 and also Johnson et al., 2014). The borehole electrodes are arranged to optimize imaging of the time-lapse distribution of changes in bulk electrical conductivity caused by the injection, migration, and reaction of the biostimulant.

Figure 15 – Configuration of ER boreholes and electrodes at the former Brandywine Defense Reutilization and Marketing Office site, Maryland, USA. Red lines (E1 through E7) represent borehole electrode arrays, blue dots represent electrodes (15 per well), and black lines (I1 and I2) show bioamendment injection locations.

Data Collection and Experimental Design

An MPT–ER¹ instrument was used to collect 10,939 dipole-dipole measurements in a variety of cross-well and in-well configurations, including full reciprocal measurements (i.e., 21,878 measurements total). The MPT-ER system uses an external voltage source for the current injection electrodes, which was set at 150 V for this survey. Current source and potential measurement electrodes were chosen according to a user-supplied measurement schedule file and switching between electrodes was completed through a control unit and a series of multiplexers attached to the electrode cables. Baseline data sets were collected before amendment injections and equivalent data sets were collected twice daily for approximately three years after injections. The amendment consisted of a proprietary lactate-based bio-stimulant with fluid conductivity significantly greater than groundwater conductivity. Two injections occurred using direct-push methods (i.e., a pipe with an opening at the end is pushed into the subsurface), one at location I1 and one at location I2 in Figure 15. For each injection, the direct-push pipe was advanced to approximately 12 m below ground surface, where approximately 113 L of amendment was injected. The pipe was then retracted by approximately 0.3 m and another 113 L was injected, repeating to approximately 1 m below ground surface. ERT data were not collected during injections but commenced shortly after the injections were completed. Importantly, the amendment formula was modified between injections such that the amendment injected into I1 had greater fluid conductivity than the amendment injected into I2.

Data Quality Control and Assurance

As previously mentioned, data quality is a function of injected current, which depends on the voltage applied to the current electrodes, electrode contact resistance, electrical conductivity of the medium, and distance between electrodes. Data quality is also a function of random and systematic error. As outlined in Section 3, random error can result from a number of factors such as: (1) high geometric factors, which result in a poor signal-to-noise ratio; (2) high contact resistance, which occurs when electrodes are in poor contact with the formation; (3) insufficient (i.e., low relative to what is needed) current injection; and (4) random ambient electrical noise. As discussed in Sections 3.2 and 3.3, random errors can be quantified from stacked and/or reciprocal measurements. It is important to eliminate (or appropriately weight) data with large random errors prior to inversion.

Histograms and summary statistics of the pre-filter reciprocity and apparent resistivity distributions for the Brandywine data are shown in Figure 16. We eliminated data with stacking errors larger than 3 percent, reciprocity larger than 10 percent, or applied currents less than 10 mA. We also eliminated data with apparent resistivities less than zero or greater than 7000 ohm-m, which would well exceed the resistivity expected for saturated medium sands at the site. Apparent resistivities above 7000 ohm-m (71 measurements) make up the tail of the histogram and were assumed to be outliers in this case (Figure 16). This data filtering was performed outside of the inversion software. The editing filters for current, reciprocity, and apparent resistivity removed 90, 58, and 751 measurements, respectively, leaving 10,040 of 10,939 measurements for a total data reduction of approximately 9 percent. Note that this type of data filtering is not required; it is possible to invert datasets with all collected data, especially if weighting each measurement individually by its data error. However, it can be useful to remove particularly poor data—the criteria for which will change with every field system such that the values here should not be assumed to be appropriate at all sites—as there is no need for the model to work to fit them.

Figure 16 – a) Reciprocal error histogram showing 10 percent data rejection limits and associated statistics. A total of 58 measurements were rejected due to excessive reciprocal error; b) Apparent resistivity histogram and rejection limits. A total of 1316 measurements were rejected due to excessive or negative apparent resistivity.

As discussed in Section 3.3, reciprocal measurements commonly are considered to provide a meaningful quantification of random errors, incorporating instrument error, error arising from nearby anthropogenic current sources, electrical storms, and/or random physical effects in violation of Equation 1 (or Equation 5). Here, however, reciprocal errors are small and have a standard deviation of 1.8 percent (Figure 16a). As shown in the next section, the reciprocal error does not represent the total error (both random and systematic as discussed in Section 4.3) and is inappropriate for use as a stopping criterion in the inversions.

Effect of Data Weighting and Regularization Weighting on Pre-Injection Inversions

Inversions were performed using a 3-D finite-element ER modeling and inversion code (E4D) developed at Pacific Northwest National Laboratory. E4D was implemented in parallel to facilitate inversion of large 3D electrical conductivity data sets. However, the Brandywine DRMO data presented here could easily be accommodated by commercially available serial inversion codes on a standard workstation. The principles and effects of data and model weighting (i.e., data noise estimation and regularization) demonstrated here are also valid for commercially available codes.

We stress that for most ER studies, the choices of data weighting and regularization will strongly affect results, and it is critical to experiment with different approaches to help distinguish artifacts from geologic features. Features common to a suite of tomograms inverted from the same dataset with different settings should inspire some degree of confidence for interpretation, whereas features that result only for a specific inversion setting or regularization are more likely to be artifacts. To demonstrate this concept, we inverted the Brandywine data under a variety of choices for data error and regularization. For the base case, we used Equation 11 to estimate data noise with a = 0.15 and b = 0.10 ohms. This choice was the outcome of trial and error using several different noise estimates, including direct reciprocal errors, and comparing inversion results with one another and with the general geology obtained from core samples. For the base-case regularization, we used an isotropic nearest-neighbor smoothing constraint, whereby neighboring elements are encouraged to be equivalent in bulk conductivity regardless of their orientation with respect to one another. The inversion started with a large ε value, and ε was decreased only when the χ² between iterations decreased by less than 5 percent. This constraint and ε ‘cooling’ procedure was chosen to ensure that any heterogeneity introduced into the inverse solution was required to fit the data (i.e., to ensure a parsimonious, or simple, solution).

To demonstrate the effects of stopping criteria on the inversion results, we show examples of inversions with different data misfits in Figure 17 and Figure 18. Note that for most inversion packages, the stopping criteria are based on some measure of data misfit as discussed previously (e.g., RMS or χ²). Because the stacking and reciprocal errors do not represent the true error (i.e., the error estimated by Equation 11), determination of the appropriate misfit is somewhat subjective. Thus, judgment and prior information are critical for determining an appropriate data fit, as indicated by the amount and distribution of heterogeneity in the inverse solution. For instance, in this case, core samples taken during installation of the ER boreholes indicated an upper layer of fine sands and silts, an intermediate zone of coarse sands and gravels constituting the aquifer, and a lower fine-grained confining unit. The approximate contact depths correspond well to electrical conductivity changes in the ‘appropriately’ fit tomogram in Figure 17b, suggesting an appropriate data weighting and stopping criteria for this inversion.

Figure 17 – A demonstration of data overfitting and underfitting. Here the data noise is estimated by Equation 11 using a = 0.15 and b = 0.1 ohm. a) Excessively smooth inverse model due to data underfitting. b) An inverse model with appropriately fit data. c) An inverse model with excessive heterogeneity due to data overfitting.

Figure 18 – L-curve plot for the Brandywine inversion. Each point and label on the plot represent an iteration of the inversion and the corresponding regularization weighting value (ε) used for that iteration. As ε decreases, more heterogeneity is introduced into the solution to decrease the χ² value, resulting in an increase in the regularization misfit. Representative visualizations of the inverse solution at selected iterations are shown to illustrate model complexity as the inversion progresses.

Figure 18 shows the L-curve representation of the inverse solution. Each point on the L-curve represents an iteration of the inversion, with the regularization term of Equation 9a and b on the horizontal axis and the χ² value (data misfit) on the vertical axis. Each point on the L-curve represents a single iteration of the inverse solution. At the homogeneous starting model, the χ² value is high, but the regularization misfit is zero. As the inversion progresses, heterogeneity is introduced into the solution to decrease the χ² value, causing an increase in the regularization misfit. When the decrease in the χ² value is small between iterations (as noted by two points close together on the L-curve), the regularization weighting (i.e., the ε value) is decreased by 50 percent, enabling the inversion to include more heterogeneity and decrease χ² at the next iteration. At the tail end of the L-curve, each decrease in ε results in a large increase in heterogeneity, but only a small decrease in the χ² value. At this point, it is assumed that the appropriate stopping point has been exceeded and the inversion is fitting noise, thus introducing artifacts into the solution. Given the data error model used in this example, a normalized χ² value of 1.08 was reached at ε = 6.25, resulting in the solution shown in Figure 17b.

In addition to stopping criteria, the data weighting scheme can have profound effects on the inversion results. Figure 19 shows two inversions with comparable data misfits [in terms of root-mean-square (RMS) error] but different data weighting schemes. Figure 19a shows results for the same data weighting scheme and regularization as the inversions in Figure 17b. Figure 19b uses data standard deviations set to a constant value of 0.1 ohm. Other than the data weighting, all inversion parameters for tomograms in Figures 19A and 19B are equal, and the inversion is allowed to iterate to an RMS error of ~0.37. The high χ² value in Figure 19b suggests that the constant standard deviation of 0.1 ohm underestimates the true error. That is, fitting the data to a χ² value of 1 in this case would overfit the data, cause artifacts, and thus yield less meaningful inversion results. However, the similarity between the tomograms for the same RMS error value suggests that the constant error estimate is useful as a relative data weight in this case.

Figure 19 – Effects of data weighting on the inverse solution. a) Data standard deviations are estimated with Equation 11 using a = 0.15 and b = 0.10. b) Data standard deviations are set to a constant value of 0.10 ohm. Each inversion was fit to approximately the same RMS value (Equation 13b). The χ² large value (Equation 13a) of inversion B suggests the data errors used for B underestimate the actual field-based errors. The differences between the inversions in a) and b) are due to the different emphasis placed on different data by the data weighting.

Inverse results also can be affected significantly by the regularization scheme, particularly in cases where the data are noisy and/or the estimated parameters are not well constrained by the data. That is, the effects of regularization are dependent on the ability of the data to resolve the electrical conductivity in all parts of the model, with resolving power varying spatially (e.g., Figure 4). Regularization will have a stronger influence on parts of the model less resolved by the data. To demonstrate the effects of regularization, we show three inversions with different regularization schemes in Figure 20. Although there are significant differences, the major features of the inversion are similar, suggesting these features are relatively well resolved by the data. It is good practice to invert data with several different regularization schemes, as this gives insight into what features are constrained by the data and what features may be controlled by regularization. Determination of which inversion features are constrained by the data and which are likely regularization-driven can also be made through model appraisal. Most model appraisal approaches require computation of the resolution matrix (Equation 15), which is only feasible for relatively small inverse problems. However, an indication of resolution is given by the sensitivity vector S (Equation 17), which shows the overall sensitivity of each region of the model to the data. The sensitivity vector S for the inverse solution in Figure 17b is shown in Figure 21. Note that the sensitivity is greater near the electrodes and in lower electrical conductivity regions. That is, a change in electrical conductivity in the resistive part of the model will have a greater influence on the data than a change in electrical conductivity in the conductive part.

Figure 20 – a) Inverse solution with isotropic nearest-neighbor smoothing constraints. b) Inversion with anisotropic nearest-neighbor smoothing constraints, encouraging homogeneity in the horizontal plane (i.e., layered structure). c) Hybrid regularization that promotes blocky structure at sharp conductivity boundaries. Data were weighted equivalent to the weighting used in Figure 14 and fit to a normalized χ² of ~1.0 and an RMS value of ~0.30 in each case.

Figure 21 – Cumulative squared sensitivity distribution (S in Equation 17) for the inverse solution in Figure 17b. Larger sensitivities are found near electrodes and in lower electrical conductivity regions.

Time-Lapse Inversions

For the time-lapse inversions, the data weighting and convergence criteria for each time-lapse data set were equivalent to that shown in Figure 17b. However, the regularization constraints were modified to encourage temporal changes in bulk conductivity with respect to the baseline inversion (Figure 17b) to vary smoothly in space. To implement those constraints, Equation 18a was added to the regularization constraint matrix, with one equation for each neighboring pair of elements in the computational mesh.

[latex]\displaystyle \epsilon (m_{i,t}-m_{i,ref})=\epsilon (m_{j,t}-m_{j,ref})[/latex]

(18a)

Here, m_i,t is the log conductivity of computational mesh element i at time t, m_j,t is a neighboring element, m_i,ref represents the pre-injection, baseline log conductivity of element i, and ϵ represents the regularization weighting parameter. Noting that:

[latex]\displaystyle m_{i,t}=m_{i,t-1}+\Delta m_{i,t}[/latex]

(18b)

(and similarly for m_j,t) Equation 18a can be re-arranged and written as Equation 18c.

[latex]\displaystyle \epsilon (\Delta m_{i,t}-\Delta m_{j,t})=\epsilon (m_{i,ref}-m_{j,ref}+m_{i,t-1}-m_{j,t-1})[/latex]

(18c)

The left- and right-hand side of Equation 18c are incorporated into the matrix equations in the left- and right-hand sides of the inversion formulation in Equation 10a, respectively.

Figure 22 shows four representative examples of the time-lapse inversion results spanning approximately one year. The elevated fluid conductivity of the amendment injected at I1 compared to I2 is immediately evident on day 5, as the increase in bulk conductivity at I1 is greater than at I2. Over time, the amendment plume spreads laterally, which is important for effective coverage and treatment. After day 93, the plume sinks and spreads over the lower flow-bounding unit, likely due to density-driven flow. By day 371, the plume is diluted relative to day 93, as expressed by a decrease in bulk conductivity.

Figure 22 – Example time-lapse inversions showing the changes in bulk conductivity caused by the presence of amendment. After injection, the amendment experienced density-driven downward flow and lateral flow over the lower bounding unit contact at approximately 10 m depth (Figure 17b). The progression from day 93 to day 371 shows the plume sinking and diluting.

¹Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.