Stratigraphic Identification with Airborne Electromagnetic Methods at the Hanford Site, Washington

Stratigraphic units can influence the fate and transport of subsurface contaminants within groundwater. Units having coarse-grained sediments act as preferential flow pathways, and therefore can accelerate the transport of contaminants to reach human and ecological receptors. At legacy waste sites, detailed knowledge of subsurface stratigraphy can be used for effective monitoring and remediation planning to help minimize risk to human health and the environment. Airborne electromagnetic (AEM) 10 methods can non-invasively provide information on kilometer-scale or larger subsurface stratigraphic features and fill informational gaps in directly sampled data from sparsely located boreholes. In this paper, we present inversion results of a 412 line-km frequency-domain AEM survey to delineate subsurface stratigraphic features at the Hanford Site, located in southeastern Washington State. The inversion was performed using a massively parallel 3D electromagnetic modeling and inversion code, where the modeling 15 is based on solving frequency-domain Maxwell’s equations using an unstructured-mesh finite-element method and the inversion employs a Gauss-Newton optimization scheme. The results are compared to an underlying geologic framework model (GFM), built by interpolating contact depths of stratigraphic units interpreted from site borehole datasets. In areas with good borehole coverage, the inversion results show a good match with the GFM to a depth of about 60 m. Outside of these areas, the inversion results exhibit 20 inconsistencies from the assumptions made to create the GFM, demonstrating that the AEM survey results can be used to improve the understanding of the geological conceptual model.


Introduction
Subsurface stratigraphy is often mapped by interpolating geologic contact depths from borehole datasets.
Near boreholes, geologic contacts are presumably more accurate than estimates at greater distances.
Interpolation of these depths may not be reliable at places where boreholes are sparsely located, or where 30 subsurface heterogeneity exists at length scales smaller than the distance between boreholes. Improved methods are needed to provide accurate spatial extent of subsurface stratigraphic units. Surface and/or airborne geophysical methods -e.g., electrical resistivity tomography (ERT), active source electromagnetic (EM) methods, seismic methods, ground penetrating radar, and magnetic methods -are obvious choices that can provide high-resolution images of the subsurface covering large areas. For 35 kilometer-scale and larger investigations, airborne EM (AEM) methods are particularly cost-effective ways to obtain subsurface coverage, with minimum disturbance, in a relatively short period of time.
There are primarily two types of AEM methods: frequency-domain and time-domain -detailed capabilities and differences of each method are presented in Siemon et al. (2009a), Steuer et al. (2009), and Thomson et al. (2007). In particular, the frequency-domain AEM method is more suitable for high-40 resolution investigations at shallower depths, while the time-domain AEM method is suitable for very large coverage areas and deeper exploration depths (Siemon et al., 2009a). Both AEM methods have been used throughout the world for mapping of geoelectrical structures in the near surface, at depths down to a few hundred meters, for the exploration of minerals (Christensen et al., 2009;Smith, 2010;Sykes et al., 2006), groundwater (Sattel and Kgotlhang, 2004;Schamper et al., 2013;Siemon et al., 2011), and oil and gas (Cox 45 et al., 2012;Pfaffhuber et al., 2009;Smith, 2010;Smith et al., 2008;Walker and Rudd, 2008); and for environmental, hazard and engineering investigations (Finn et al., 2007;Goebel et al., 2019;Minsley et al., 2012;Pfaffhuber et al., 2017). These applications utilize the fact that the interpreted electrical resistivity (or conductivity) anomalies can directly be correlated with subsurface properties. For example, in sedimentary rocks, the resistivity tends to decrease as the grain size decreases and as clay content increases; 50 unsaturated sequences will have a higher resistivity value than their saturated counterparts; hydrocarbon saturation in pore space tends to increase resistivity compared to water saturation (Palacky, 1987).
As AEM surveys generate a large amount of data consisting of tens of thousands of sounding locations, data are normally interpreted using 1D analyses, e.g., by applying apparent resistivity transforms (Huang and Fraser, 1996;Sengpiel, 1988), conductivity-depth imaging (Huang and Rudd, 2008;Macnae et al., 55 1991), 1D inversions (Chen and Raiche, 1998;Farquharson et al., 2003), or laterally/spatially constrained 1D inversion (Auken and Christiansen, 2004;Siemon et al., 2009b;Viezzoli et al., 2008). Despite widespread acceptance and the argument that spatially constrained 1D inversion can provide good results for geological settings where each observation can be simulated with 1D forward modeling Viezzoli et al., 2010), it has been demonstrated that such 1D analyses are usually not adequate when 60 2D or 3D geological complexity exists (e.g., see Cox et al., 2010;Ellis, 1998;Raiche et al., 2001;Wilson et al., 2006). Several 2D and 3D modeling and inversion analyses have also been developed (e.g., by Cox et al., 2010;Heagy et al., 2017;Noh et al., 2018;Sasaki, 2001;Xi and Li, 2016;Yang et al., 2014); however, their acceptance has been limited because the number of observations used in the inversion was restricted due to computational constraints, or because the inversion codes were not publicly available.

65
In 2008, frequency-domain AEM data were acquired in the eastern part of the Hanford Site to delineate highly transmissive zones that may serve as preferential flow pathways for subsurface contaminants. Initial interpretation of the data was based on a 1D analysis by transforming the flight-line data into depth slices of electrical resistivity by computing apparent resistivity and apparent depth (Huang and Fraser, 1996).
This transformation showed a weak correlation of resistivity with the underlying stratigraphy derived from 70 nearby boreholes, which is likely the result of representing the subsurface from interpolation of simplified 1D analysis as opposed to a 2D or 3D analysis. In 2010, two ground-based geophysical surveys, ERT and frequency-domain ground-based EM, were conducted along a test profile adjacent to a portion of the AEM survey flight lines, with the main objective of validating the AEM data interpretation (CHPRC, 2010a(CHPRC, , 2010b. Both ERT and ground-based EM data were inverted using 2D inversion schemes and the obtained 75 resistivity models were fairly consistent with the stratigraphic units interpreted from six nearby borehole datasets down to 30-40 m depth. However, these inversion results showed poor agreement with the resistivity slices computed from the AEM data. In this paper, we present the first inversion results of the AEM data at the Hanford Site obtained using a 3D EM modeling and inversion code. This code has recently been developed as a part of a capability extension 80 to E4D, an open-source software widely used for 3D/4D ERT data modeling and inversion (Johnson et al., 2010(Johnson et al., , 2017Johnson and Wellman, 2015). The forward modeling is based on solving frequency-domain Maxwell's equations using an unstructured-mesh finite-element method, and the inversion uses a Gauss-Newton optimization scheme. The implementation is massively parallelized, which allows for forward and inverse simulations to be efficiently distributed on thousands or even more processing cores of a 85 supercomputer for large-scale problems. For survey configurations similar to the collected AEM data, a synthetic inversion study was first executed to determine if a full 3D inversion could recover the subsurface resistivity between flight paths. Inversions were then performed along individual flight lines because the synthetic study demonstrated that the flight-line spacing was too large for a full 3D inversion. The inversion results were subsequently compared to an existing geologic conceptual model developed by interpolating 90 contact depths of stratigraphic units interpreted from site borehole datasets.
The remainder of this paper is organized as follows: we first present a brief description of our study area, the Hanford Site. We then describe our methodology with details on the AEM system and data acquisition, 3D EM modeling and inversion code, and relationship of electrical resistivity with different stratigraphic units at the Hanford Site. We thereafter present inversion results of the AEM data, followed by their 95 interpretation. Finally, we draw concluding remarks on our work.

Study area
The Hanford Site is a part of the U.S. Department of Energy (DOE) nuclear weapons complex. It is located  (Gephart, 2003;Lichtenstein, 2004). According to published historical reports on the Hanford Site, 67 single-shell tanks have, or are suspected to have, collectively leaked over 1 110 million gallons of highly radioactive wastes into the vadose zone (Gephart and Lundgren, 1998;Gephart, 2003;Lichtenstein, 2004).
The hazardous wastes presented a risk of contaminating groundwater and possibly the nearby Columbia River through subsurface flow pathways. A recent groundwater monitoring report at the Hanford Site reported that about 168 km 2 of groundwater is contaminated above regulatory standards (DOE, 2019).

115
Groundwater across the Hanford Site typically flows from west to east and discharges along the Columbia River. Stratigraphic units and paleochannels, filled with coarse-grained sediments, act as preferential flow pathways to contaminants within groundwater, and consequently can shorten travel times to reach human and ecological receptors. It has been estimated that the contaminant travel time from the Central Plateau to the Columbia River can vary from a few decades for the most mobile contaminants to a century or more 120 for the less mobile ones within the groundwater (Gephart, 2003). For effective monitoring and remediation planning, it is therefore critical to have an accurate understanding of the stratigraphic units and possible flow pathways that can influence the subsurface contaminant transport across the Hanford Site.
The generalized stratigraphic units of the Hanford Site consist of, in descending order, the Hanford formation, the Ringold Formation, and the Columbia River Basalt Group (Fig. 2). This simplified layer-125 cake depositional model has been complicated by the depositional process, and subsequent removal, of sedimentary units. During the Ringold Formation depositional cycle, actions of ancestral rivers resulted in deposition of intercalated layers of indurated to semi-indurated sediment including clay, silt, fine-to coarsegrained sand, and granule-to-cobble gravel. Within most of the Ringold Formation, there is a thick confining layer of lower mud, called the Ringold Mud unit, composed of lacustrine silts and clay, paleosol, 130 and fluvial overbank. After the Ringold Formation deposition, regional incision occurred for an extended period, followed by soil formation and deposition of windblown sediments. These sedimentary deposits are termed the Cold Creek unit. This was followed by the deposition of coarse gravel-to-boulder sediments of the Hanford formation due to cataclysmic floods during the last ice age. The floods also incised channels of different depth into the previously deposited Ringold and Cold Creek sediments, and in some cases,

135
probably the basaltic bedrock. These scoured areas and paleochannels were later filled by poorly sorted Hanford sands and gravels. The above descriptions are gathered from DOE (2011); however, for comprehensive details on the stratigraphy at the Hanford Site, interested readers are further referred to DOE (2002), Lindsey (1995), and Lindsey et al. (1994).
Across the Hanford Site, these stratigraphic units have continuously been mapped to build a 3D geologic 140 framework model (GFM) by interpolating depths of different geologic contacts interpreted from approximately 1500 boreholes (CHPRC, 2018). Borehole datasets may include drillers' descriptions, geologists' descriptions, geophysical logs, grain-size analyses, and/or sediment photographs, and were also used to establish the presence of ancestral fluvial channels or paleochannels across the site (Bjornstad et al., 2010;Reidel and Chamness, 2007). The 3D GFM is expected to be accurate in areas with higher 145 borehole densities; however, it may have more uncertainty in areas with sparsely located boreholes or areas where subsurface heterogeneity exists at length scales smaller than the borehole separation.

AEM system and data acquisition
system was used to collect AEM data. The objective of this AEM survey was to delineate stratigraphic/lithologic changes in the upper layers to a depth of approximately 50 m. The RESOLVE ® system is based on the frequency-domain AEM method and measures the secondary EM field induced in the ground by a transmitter loop. It operates with five pairs of horizontal co-planar (HCP) loop-loop configurations at 400, 1800, 8200, 40000, and 140000 Hz nominal frequencies, and one pair of 155 vertical co-axial (VCX) loop-loop configuration at 3300 Hz. The transmitter-receiver separation is 7.9 m for the HCP coil sets and 9.0 m for the VCX coil set. The transmitter and receiver coil sets are mounted in a "bird", which is towed beneath an aircraft or helicopter during data acquisition.
The blue enclosed area, adjacent to the Columbia River, in Fig. 1 shows the AEM survey area at the Hanford Site 600 Area. The RESOLVE ® system was towed beneath a helicopter and flown at an average speed of 160 133 km/h and a nominal ground clearance of 30 m over the survey area between June 29 and July 1, 2008.
In some portions of the survey area, e.g., above power lines, the ground clearance was much higher than the nominal ground clearance due to safety reasons. The actual operational frequencies were 385, 1792,

AEM data modeling and inversion
The forward modeling presented in this paper is performed using a 3D unstructured-mesh finite-element forward modeling software. Assuming a temporal-dependence of e !"#$ with the angular frequency " and 185 i = √−1, a quasi-static vector Helmholtz equation for the electric field can be derived from Maxwell's equations. Following Jaysaval et al. (2014Jaysaval et al. ( , 2015Jaysaval et al. ( , 2016 where e is the intensity of the electric field; -is the electric field source; and μ and σ are, respectively, the magnetic permeability and electrical conductivity of the medium. The electrical conductivity is inverse of the electrical resistivity /.

190
To get accurate results for a frequency-domain AEM system, we solve Eq. (1) by decomposing total electric field ) into primary ) % and secondary ) & fields (Newman and Alumbaugh, 1995) as and where 2 % is the background conductivity of the medium. Using Eqs.
(1), we get The edge based finite element discretization of Eq. (4) on a tetrahedral mesh is assembled into a system of linear equations (Jin, 2002): where 4 is the system matrix, 5 is the vector of unknown electric field, and 6 is the source vector resulting from the right-hand side of Eq. (4). This equation is solved to compute the electric field followed by using Faraday's law to compute the magnetic field. Accuracy of the developed forward modeling software was 200 benchmarked against analytical/semi-analytical solutions for layered earth and Jaysaval et al. (2014Jaysaval et al. ( , 2015Jaysaval et al. ( , 2016 for 3D earth models. The results agreed well, within acceptable error ranges, depending on fine-or coarse-meshing of benchmarking models. In the inversion, we have implemented a Gauss-Newton optimization scheme following Johnson et al. (2010). The objective of the inversion is to find a conductivity model 7 (in our case it represents 205 logarithmic of the conductivity) such that it minimizes a cost function 8: where 8 ' is the data cost function measuring the misfit between the observed : )*& and forward modeled : &+, data; 8 ( is the corresponding misfit measure between 7 and constraints placed upon it; and 9 is a regularization parameter controlling the contribution of 8 ( to 8 compared to 8 ' . We choose to minimize 8 using the least-square method, and hence the data and model cost functions are, respectively, and Here, 7 -./ is the logarithmic of a reference conductivity model, < ' is the data-weighing matrix, and < ( is the regularization matrix.
We set ∂8 ∂ ⁄ 7 = 0 to minimize 8 using the least-square method. This results in the following normal equation: where I is the Jacobian or sensitivity matrix that represents partial derivatives of the synthetic data with respect to the model parameters and δ7 is a model update vector. The superscripts T and H denote, respectively, the transpose and conjugate transpose operators. The inversion is performed iteratively: at each iteration with a current model 7 2 , we seek to obtain δ7 2 by solving normal Eq. (9) such that a new model 220 7 234 = 7 2 + N δ7 2 decreases the cost function 8. Here, N is the step length. This iterative process continues until the data misfit reaches below a pre-defined tolerance level.
The Jacobian matrix is computed using the adjoint-state method (McGillivray et al., 1994), and the normal equation (Eq. 9) is solved using a conjugate-gradient least-squares method following Johnson et al. (2010).
Further implementation details of the Jacobian matrix and the solver are beyond the scope of the current 225 paper. However, we do emphasize that they are massively parallelized and usually computed using several thousands of processing cores for large-scale 3D problems.

Electrical resistivity and lithology
The application of any geophysical EM or electrical method is primarily based on the relationship of electrical resistivity with lithology, porosity, and/or pore fluids (Palacky, 1987). The knowledge of such a 230 relationship will guide us to interpret resistivity models obtained by inverting the collected AEM data at the Hanford Site 600 Area.
In sedimentary rocks, there are two dominant mechanisms that contribute to electrical resistivity: ionic conduction and surface conduction (Knight and Endres, 2005). The ionic conduction occurs through pore fluids and its efficiency depends on porosity, pore-fluid saturation, and pore-fluid type, e.g., in sedimentary 235 rocks with water in pore space, the resistivity tends to decrease as water saturation increases. As a result, fully saturated sedimentary sequences at the Hanford Site will have a lower resistivity value than their lesssaturated counterparts. Ionic conduction is also affected by the grain size: the coarser the grain size, the higher will be electrical resistivity because the bulk of the volume is dominated by insulating minerals.
Therefore, gravels and sands of the Hanford formation are expected to have higher resistivity values than 240 adjacent finer-grained sediments of the Ringold Formation.
The second mechanism, the surface conduction, occurs due to the presence of a high concentration of ions associated with the electrical double layer at the solid/water interface (Knight and Endres, 2005). The occurrence of surface conduction causes a decrease in electrical resistivity and its efficiency depends on the surface area: the higher the surface area, the higher the surface conduction, which is inversely related 245 to the electrical resistivity. As clays and finer grains have higher surface areas, the resistivity tends to decrease as the clay content increases and as the grain size decreases. Consequently, at the Hanford Site, the finer-grained sediments of the Ringold Formation are expected to have lower resistivity values than coarse-grained sediments. For example, sediments in the Ringold Mud unit -zones of increased silt/clay contents within the Ringold Formation -are expected to have low resistivity values.

250
In igneous rocks, the resistivity depends on their silica content: the more silica the rock matrix contains, the more resistive the rock becomes (Palacky, 1987). Although basalt has relatively less silica content compared to other igneous rocks, e.g., felsic igneous rocks, its resistivity is usually much higher than common sedimentary rocks and can vary in the order of 1000 Ωm. The simulated AEM data were inverted using the previously described inversion algorithm using a starting model having half-space resistivity of 500 Ωm. The inversion iterations were stopped after the data misfit reduced to 1% of the initial data misfit value in all the three cases. Figs (Fig. 4c), the inversion recovered the true resistivity, but produced some artifacts, shown as intervals of highly resistive features in the top layer.

Inversion of AEM data at the Hanford Site 600 Area
In line. The inversions were allowed to converge when the data misfit reduced to about 5% of the initial data misfit value. Plots comparing the in-phase and quadrature components of the observed and predicted data 295 for the inverted resistivity models are shown in Fig. 5 for all five frequencies. These plots suggest that the predicted data are qualitatively agreeing well with the observed ones.
We now show the inverted resistivity model along flight-line L10300 (see Fig. 3 for the line orientation) in Fig. 6. The model is shown down to 100 m depth. However, the result may be reliable only up to about 60 m depth, which is the effective penetration depth limit of the RESOLVE ® AEM system reported by the 300 operator (Fugro). Compared to other flight lines, this line has the most boreholes nearby with good coverage; it has three boreholes, B-15, B-16, and B-17 (shown by green dots near line L10300 in Fig. 3 and by black vertical lines in Fig. 6), which helps to validate the inversion result. Three geologic contacts, the top of the Ringold E unit, the Ringold Mud unit, and the Columbia River Basalt Group (hereinafter the basalt), extracted from the GFM are overlain on this resistivity model and shown, respectively, by black, 305 magenta, and orange curves in Fig. 6. These three contacts were also interpreted from boreholes B-15, B-16, and B-17, and are marked, respectively, by black, magenta, and orange cross marks. We notice a minor mismatch between some of the geologic contacts derived from the boreholes and GFM. This is because the boreholes may not be located exactly below the flight line. Two power lines passing this cross-section are marked by P-1 and P-2. These power lines yielded wide swaths of unreliable AEM data, and therefore   Table   1) for lines L10370-L10310 correlate fairly well with the geological contacts derived from GFM and boreholes. Note that flight-lines L10300-L10370 have reasonably good coverage of boreholes; therefore, the interpolated GFM below these flight lines is expected to be less uncertain.
Figs. 8a-8i present the inverted resistivity models, respectively, along flight-lines L10150-L10070. As 345 done previously, the geologic contacts and power lines are shown on these inverted models. Below these flight lines, there is no borehole coverage; therefore, the GFM may be more uncertain as compared to those shown in Figs. 6 and 7. The Hanford-Ringold E contact derived from the GFM correlates well with the inverted resistivity model (see Table 1) for the easternmost flight-line L10150. However, as we move westward to line L10070 from line L10150, this contact appears to go deeper in the inverted resistivity 350 models than the one derived from the GFM. Nonetheless, the Ringold E-Ringold Mud contact derived from the GFM seems to agree reasonably well with the inverted resistivity models.

360
Ringold E contact appears to be deeper relative to the borehole B-12 derived contact (see Fig. 9i).
Nevertheless, B-12 is located very close to power line P-1 and, as mentioned previously, the inverted resistivity values should not be considered for the interpretation. Generally, for these flight lines, the geologic contacts derived from GFM show poor agreements with the ones interpreted from the inverted

375
We have successfully inverted 412 line-km of frequency-domain AEM data at the Hanford Site with a massively parallel 3D EM modeling and inversion code. The inversion results delineated subsurface stratigraphic units at the Hanford Site. In areas with good coverage of boreholes, the results showed a good match with an existing GFM built by interpolating geologic contact depths derived from these boreholes.
Outside these areas, where the framework model is more uncertain, the inversion results provided additional 380 information on geologic contact locations for the stratigraphic units. This knowledge can further be incorporated into the GFM to reduce its uncertainty and provide a technical basis for monitoring and remediation activities.
Author contributions. JLR and PJ planned the study. PJ implemented the modeling and inversion code, performed the inversion, analyzed the results, and made the figures. JLR contributed to analyze the results

385
and to make figures for the location map. TCJ was supervising the study. PJ wrote the manuscript, and all co-authors read the manuscript and helped with editing.
Competing interests. The authors declare that they have no conflict of interest.