3D inversion of airborne electromagnetic data

Time-domain airborne surveys gather hundreds of thousands of multichannel, multicomponent samples. The volume of data and other complications have made 1D inversions and transforms the only viable method to interpret these data, in spite of their limitations. We have developed a practical methodology to perform full 3D inversions of entire time- or frequency-domain airborne electromagnetic (AEM) surveys. Our methodology is based on the concept of a moving footprint that reduces the computation requirements by several orders of magnitude. The 3D AEM responses and sensitivities are computed using a frequency-domain total field integral equation technique. For time-domain AEM responses and sensitivities, the frequency-domain responses and sensitivities are transformed to the time domain via a cosine transform and convolution with the system waveform. We demonstrate the efficiency of our methodology with a model study relevant to the Abitibi greenstone belt and a case study from the Reid-Mahaffy test site in Ontario, Canada, which provided an excellent practical opportunity to compare 3D inversions for different AEM systems. In particular, we compared 3D inversions of VTEM-35 (time-domain helicopter), MEGATEM II (time-domain fixed-wing), and DIGHEM (frequency-domain helicopter) data. Our comparison showed that each system is able to image the conductive overburden and to varying degrees, detect and delineate the bedrock conductors, and, as expected, that the DIGHEM system best resolved the conductive overburden, whereas the time-domain systems most clearly delineated the bedrock conductors. Our comparisons of the helicopter and fixed-wing time-domain systems revealed that the often-cited disadvantages of a fixed-wing system (i.e., response asymmetry) are not inherent in the system, but rather reflect a limitation of the 1D interpretation methods used to date.

3D inversion of airborne electromagnetic data Leif H. Cox1, Glenn A. Wilson1, and Michael S. Zhdanov2 ABSTRACT Time-domain airborne surveys gather hundreds of thousands of multichannel, multicomponent samples. The volume of data and other complications have made 1D inversions and trans-forms the only viable method to interpret these data, in spite of their limitations. We have developed a practical methodology to perform full 3D inversions of entire time- or frequency-do-main airborne electromagnetic (AEM) surveys. Our methodol-ogy is based on the concept of a moving footprint that reduces the computation requirements by several orders of magnitude. The 3D AEM responses and sensitivities are computed using a frequency-domain total field integral equation technique. For time-domain AEM responses and sensitivities, the frequency-domain responses and sensitivities are transformed to the time domain via a cosine transform and convolution with the system waveform. We demonstrate the efficiency of our methodology with a model study relevant to the Abitibi greenstone belt and a case study from the Reid-Mahaffy test site in Ontario, Canada, which provided an excellent practical opportunity to compare 3D inversions for different AEM systems. In particular, we com-pared 3D inversions of VTEM-35 (time-domain helicopter), MEGATEM II (time-domain fixed-wing), and DIGHEM (fre-quency- domain helicopter) data. Our comparison showed that each system is able to image the conductive overburden and to varying degrees, detect and delineate the bedrock conductors, and, as expected, that the DIGHEM system best resolved the conductive overburden, whereas the time-domain systems most clearly delineated the bedrock conductors. Our comparisons of the helicopter and fixed-wing time-domain systems revealed that the often-cited disadvantages of a fixed-wing system (i.e., response asymmetry) are not inherent in the system, but rather reflect a limitation of the 1D interpretation methods used to date. INTRODUCTION To improve mineral exploration success, there is industry-wide consensus on the need to increase the "discovery space" by explor-ing under cover and to greater depths. Terranes characterized by laterally varying conductive near-surface geology produce EM re-sponses which cannot be accounted for with 1D inversion. These can then produce artifacts which manifest themselves as conductive features at depth, effectively masking the more resistive basement and deeper conductors, which potentially contain economic miner-alization. The exploration industry has responded to this by devel-oping airborne electromagnetic (AEM) systems with ever-higher transmitter moments, and hardware and processing technology have improved data quality significantly. For example, high-moment AEM systems such as MEGATEM II (Smith et al., 2003) and SPECTREM 2000 (Leggatt et al., 2000) claim maximum depths of penetration up to 900 m, subject to terrain conductivity, noise, and system specifications. While AEM systems continually evolve with higher moments, improved geologic modeling through 3D inversion has the potential to increase the success of AEM-led exploration. The Abitibi green-stone belt, for example, is a world-class volcanogenic massive sul-fide (VMS) district with a statistically significant population of yet-to- be-found VMS deposits. Following the MEGATEM-led discov-ery of Perseverance deposit in the Matagami camp in 2000, no further major discoveries have been attributed to the 180; 000þ line km MEGATEM initiative of the 2000s, where interpretations were based on empirical analysis or 1D inversion of the observed AEM data. At the time, modeling by Noranda (now Xstrata) de-monstrated that the MEGATEM system could detect VMS deposits to at least 250 m depth (Smith et al., 2003). Hence, the lack of Manuscript received by the Editor 3 October 2011; revised manuscript received 9 April 2012; published online 10 July 2012. 1TechnoImaging, Salt Lake City, Utah, USA. E-mail: [e-mail]; [e-mail]. 2University of Utah, Department of Geology and Geophysics, Salt Lake City, Utah, USA; TechnoImaging, Salt Lake City, Utah, USA. E-mail: mzhdanov@ technoimaging.com. © 2012 Society of Exploration Geophysicists. All rights reserved. WB59 GEOPHYSICS, VOL. 77, NO. 4 (JULY-AUGUST 2012); P. WB59-WB69, 12 FIGS., 1 TABLE. 10.1190/GEO2011-0370.1 Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ discoveries does not represent any limitation of the MEGATEM system. Rather, it can be argued that the lack of discoveries demon-strates the limitations on how the MEGATEM data has been inter-preted to date. This suggests that potential exists for making new discoveries in the Abitibi greenstone belt by reinterpreting the his-toric MEGATEM data with 3D inversion. Similar conclusions can be drawn elsewhere in the world. That said, to influence exploration decisions in a tangible man-ner, 3D AEM inversion is required to invert entire surveys (>1000 line km) with a deposit-scale model resolution (e.g., ≈ 25 m) in sufficient time. Until recently, there was no practical 3D AEM inversion available in either industry or academia. To quote Macnae (2007), "progress toward routine 2D and 3D inver-sions of AEM data has been slow, despite significant effort." A his-tory of 3DAEM inversion is summarized in Figure 1. As a metric of inversion capability, we define the 3D AEM problem size as the product of the number of cells in the 3D model and the number of AEM stations modeled. While 3D modeling methods and their approximations have varied between authors (e.g., Ellis, 1995, 2002; Zhdanov and Tartaras, 2002; Zhang, 2003; Zhdanov and Chernyavskiy, 2004; Raiche et al., 2007; Oldenburg et al., 2010), 3D AEM inversion has been limited to relatively small kilo-cell 3D models with tens of stations. Again, to quote Macnae (2007), "while the formal inversion process is practical for targets already identified as being of interest, it is unlikely to be routine in the near future as a means of processing complete surveys." Cox and Zhdanov (2006) used a combination of the full integral equation (IE) method for modeling, and localized quasilinear (LQL) method for calculating the sensitivities. Cox and Zhdanov (2007) used the same combination of the full IE and LQL methods, yet they introduced the novel concept of a moving footprint for computing the responses and sensitivities. By doing so, Cox and Zhdanov (2007) were able to increase the AEM problem size by nearly five orders of magnitude. Using a similar moving footprint methodology, Cox et al. (2010) used the full IE method for both modeling and calculating sensitivities. This clearly represents a paradigm change in 3D AEM inversion methodology, and unlike others, it has enabled practical 3D inversion of entire surveys. In the present paper, we extend our method for a frequency-domain AEM problem, detailed in Cox et al. (2010), to time-domain AEM surveys and compare the results of 3D inversion of both frequency-domain and time-domain AEM data. We base our frequency-domain modeling on the 3D contraction integral equation method (Hursán and Zhdanov, 2002). For time-domain AEM, system responses and sensitivities are obtained by Fourier transform of the frequency-domain responses and sensitivities, and are convolved with the transmitter waveform (Raiche, 1998). This enables us to accurately model and invert data from any extant or proposed AEM system. We use a regularized conjugate gradient method for minimizing our objective functional with the smooth or focusing stabilizers (Zhdanov, 2002). Our implementation makes it practical to invert entire surveys with thousands of line-km of AEM data to megacell 3D models using a multiprocessor workstation (e.g., black diamonds in Figure 1). In the public domain, case stu-dies have been presented for VTEM (e.g., Wilson et al., 2011a; Combrinck et al., 2012), SPECTREM (e.g., Pare et al., 2012), TEMPEST (e.g., Wilson et al., 2011b), SkyTEM (Wilson et al., 2010) and RESOLVE (Cox et al., 2010) surveys. In this paper, we demonstrate the effectiveness of our 3D AEM inversion with a synthetic model study relevant to the Abitibi greenstone belt, and follow with a case study from the Reid-Mahaffy test site in On-tario, Canada, where we have inverted DIGHEM, VTEM, and MEGATEM data to 3D conductivity models that independently capture the geology as known from drilling. 3D MODELING AND INVERSION METHODOLOGY The principles of regularized 3D AEM inver-sion were outlined in Zhdanov (2009). However, at each iteration of a nonlinear inversion, compu-tation of the responses and their sensitivities for an entire AEM survey is not trivial. Cox et al. (2010) introduce a moving footprint methodol-ogy to minimize these calculations for the fre-quency- domain AEM, and we refer the reader to that paper for a detailed description of our moving footprint methodology. In what follows, we present time-domain implementation of this method. Modeling Time-domain AEM modeling can be accom-plished either by direct, time-domain solutions or by transformation of frequency-domain solu-tions. The latter offers distinct advantages; namely, that the effects of frequency-dependent conductivity (e.g., induced polarization) can be modeled, that artificial dispersion effects that Figure 1. A brief history of 3D AEM inversion from 1995 to the present, summarized in terms of AEM problem size, i.e., number of cells in the 3D model times the number of stations in the AEM survey. Cox and Zhdanov (2007) achieved nearly a five order of magnitude increase over Cox and Zhdanov (2006) by introducing a moving footprint. Black diamonds represent a subset of the 3D AEM inversions completed by the authors. WB60 Cox et al. Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ arise in direct time-domain solutions are avoided, and that the many linear system solutions required for implicit time stepping are avoided. The latter point is particularly critical for accurately mod-eling AEM systems with nonideal transmitter waveforms. As described in Cox et al. (2010), our modeling is based on the rigorous 3D contraction IE method (Hursán and Zhdanov, 2002). We do not use any Born or Born-like linear approximations in our modeling of the fields or their sensitivities. While the IE method has been detailed elsewhere, we do provide a short description of our particular implementation of it that has not been described pre-viously. In discrete form, the IE method requires the solution to the following system of equations for each sensitivity domain (foot-print) ðI − ΓE · σaÞ · Ea ¼ ΓE · σa · Eb; (1) where Ea is the vector of the anomalous a electric field, I is the identity matrix, ΓE is the matrix of volume-integrated electric Green's tensors for the background b conductivity model, and σa is a diagonal matrix of anomalous conductivities. Eb is the vector of the background electric fields. Each footprint domain is a self-contained domain with respect to the forward modeling. It is a rec-tangular anomalous domain with sides the size of the footprint, ex-tending to the bottom of the inversion domain, and embedded in a layered earth. This domain is no different than that which would be set up for forward modeling a single transmitter-receiver pair. For each footprint domain, the background conductivity model is cho-sen such that the background magnetic fields will dominate the AEM system response at both early and late times. The background conductivity model is horizontally layered so that the background electric and magnetic fields can be solved semianalytically in a stable and accurate manner. We note that for AEM systems, there is negligible error in representing the transmitter as a unit magnetic dipole rather than a polygonal loop due to the rapid lateral expan-sion of the "smoke ring" in the air prior to interaction with the earth (Raiche, 1998). Because all footprint domains are equally discre-tized, the electric Green's tensors ΓE need only be calculated once, and then translated over the 3D earth model. The magnetic fields at the receivers are then computed from H ¼ Hb þ ΓH · σa · ½Eb þ Ea; (2) where ΓH is the matrix of volume-integrated magnetic Green's ten-sors for the background b conductivity model. Note that the AEM transmitter and receiver pitch, roll, and yaw can be included via appropriate Euler rotations of the background electric and magnetic fields and the magnetic Green's functions. Equation 2 requires the total electric field in each cell, and this is computed as the sum of the background and anomalous electric fields. For models with high conductivity contrasts or very resistive hosts, the background and anomalous electric fields are of near-equal amplitude but op-posite sign. Given finite precision, their addition introduces numer-ical errors for conductivity contrasts greater than 300:1 that are artificially inflated by equation 2 (e.g., Xiong et al., 1999). By add-ing the background electric fields to both sides of equation 1 and with some algebra, we obtain the linear system ðI − ΓE · σaÞ · E ¼ Eb: (3) This means we can solve for the total electric field instead of the anomalous electric field, and still retain the same distributed source term (Hursán and Zhdanov, 2002). This type of formulation is un-ique to integral equations methods and their hybrids. Over the past decade, our IE method was carefully verified and tested by members of the University of Utah's Consortium for Elec-tromagnetic Modeling and Inversion (www.cemi.utah.edu). As an illustration of the accuracy of our IE method, we present in Figure 2 a comparison between the impulse response calculated from our 3D IE method and a 1D semianalytic method for a HELITEM time-domain helicopter system above a layered half-space. The 1D model contained three layers, with an upper 100 m thick layer of 3000 ohm-m, a middle layer 100 m thick of 1 ohm-m, and a basement of 3000 ohm-m, so a resistivity contrast of 3000:1 was simulated. As shown in Figure 2, the 3D IE results match the 1D semianalytic results very well. We should note that the anom-alous fields were very large in this example due to the high contrast with background, yet the 3D IE method produced very accurate re-sults. This confirms once again that our 3D IE method is accurate for most mining applications, and virtually all environmental and hydrocarbon applications. For calculating time-domain AEM responses, it is not a simple matter of fast Fourier transforming a frequency-domain spectrum. Using the method of Raiche et al. (2007), the frequency-domain response is computed at six frequencies per decade, logarithmically spaced from 10 mHz to 1 MHz, that correspond to the notches in a digital filter. The imaginary components are splined and extrapo-lated back to zero frequency. The time-domain step response, Bðr; tÞ; is then calculated from the cosine transform Bðr; tÞ ¼ − 2μ0 π Z ∞ 0 ℑ½Hðω; rÞ ω cosðωtÞdω; (4) which is evaluated using a digital filter. The impulse response is not computed directly due to singularity at zero time. The step response is computed to several pulse lengths, and is then folded and stacked into a single pulse length and resplined. This emulates stacking dur-ing acquisition. The system step response is then obtained by con-volving the spline of the derivative of the transmitter current ∂IðtÞ∕∂t with the step response Figure 2. Comparison of 1D semianayltic (shown by a black line) and 3D integral equation (shown by crosses) responses for a HE-LITEM time-domain helicopter system above a layered earth with 1:3000 resistivity contrast. 3D AEM inversion WB61 Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ Sðt; rÞ ¼ ∂IðtÞ ∂t Bðr; tÞ. (5) The system impulse response can be obtained from analytic differ-entiation of the spline coefficients of the system step response. The system responses are then numerically integrated over the receiver windows. Sensitivities are computed in the frequency-domain using adjoint operators (e.g., Zhdanov, 2002) and transformed to the time domain AEM system sensitivities as just described. In practice, the transmitter current waveform and primary field is measured by flying the AEM system to a high altitude. If these data are available, the primary is often subtracted from the survey data, with the result normalized to the peak value of the primary field. This primary field differs from the ideal DC primary field because it contains the aircraft's time-domain response as well. The fact that the observed and predicted primary fields differ can introduce errors in inversion. This problem is overcome by convolving the measured primary field or calibration waveform rather than the current deri-vative and normalizing to the peak value of this same waveform (Raiche, 1998). If the primary field calibration is not available, the transmitter current waveform is used. Inversion AEM surveys are finite in their spatial and frequency content and contaminated with noise, and their inversion is inherently nonuni-que and unstable, meaning it is an ill-posed problem. Regularization must be introduced so as to obtain a unique and stable solution. Regardless of the methodology used, most regularized inversion seeks to minimize the Tikhonov parametric functional, PαðσÞ PαðσÞ ¼ ϕðσÞ þ αkWmðσ − σapr Þk2 → min : (6) In the last equation, α is the regularization parameter, and ϕðσÞ is a misfit between the predicted, dpred ¼ AðσÞ; and observed, dobs; data, ϕðσÞ ¼ kWdðAðσÞ − dobs Þk2; (7) where A is the nonlinear forward operator, σ is the Nm length vector of conductivities, dobs is the Nd length vector of observed data, σapr is the Nm length vector of a priori conductivities, and k : : : k denotes the respective Euclidean norm. The Tikhonov functional is mini-mized in the same manner as Cox et al. (2010). Data and model weights are introduced to equation 6 through data and model weighting matrices, Wd and Wm, respectively, which weight the inverse problem in logarithmic space to reduce the dynamic range of the data and conductivities. Our model weights are simply the inverse of the square root of the integrated sensitivity for each cell (Zhdanov, 2002). We use different data weights for time- and fre-quency- domain AEM data. Frequency-domain AEM data can vary by orders of magnitude with frequency and component. For example, for the same fre-quency, horizontal coplanar data are typically larger than vertical coaxial data. Moreover, vertical coaxial data are noisier than hor-izontal coplanar data. It follows that for every frequency and com-ponent of frequency-domain AEM data, we can use the data weights Wd;i ¼ 1 ϵi ; (8) where ϵi is the error at each data point. This data weight can be generalized to vary with stations, though acquisition contractors ty-pically quote errors only for frequencies and components. These data weights take into account the different noise levels for the var-ious frequencies and components. Time-domain AEM data can vary by several orders of magnitude, and noise levels can be quoted as a single value for all channels, components, and stations. Hence, time-domain AEM inversion with data weights as per equation 8 will preferentially fit the large-amplitude early-time channels rather than the smaller-amplitude late-time channels. It follows that for every channel, component, and station of time-domain AEM data, we instead use the point symmetric data weights Wd;i ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 2 ðd2 obs;i þ d2 pred;i Þ s ; (9) such that the residual errors across all channels, components, and stations are weighted equally. This means that the residual errors at a station above an anomaly of interest are just as important as the residual errors somewhere removed from the region of interest, or that residual errors at the earliest times are just as important as those at the latest times. This means the misfit calculated using equation 9 is usually higher than a misfit calculated from equation 8 for the same observed and predicted data. With these data weights, the re-sidual errors are also symmetric with respect to the observed and predicted data. Moreover, this data weight bounds the misfit be-tween 0% and 200%; the maximum misfit occurring when observed and predicted data have the same amplitude but opposite sign. Using this type of data weights does require a noise floor, so data below the instruments noise level or near zero crossing are not given high importance.We specify a noise floor, below which the data are simply rejected. With a moving footprint, the Fréchet matrix can be constructed as a sparse matrix with memory and computational requirements re-duced by several orders of magnitude. The number of nonzero ele-ments in each row of the sensitivity matrix is just the number of elements within each footprint (on an order of thousands) rather than the total number of elements in the domain (in the order of hundreds of thousands or millions). As expected, helicopter fre-quency- domain systems such as RESOLVE and DIGHEM have the smallest footprints. Helicopter time-domain systems such as SkyTEM, VTEM, AEROTEM, and HELITEM have larger foot-prints. Fixed-wing time-domain systems such as TEMPEST, SPEC-TREM, GEOTEM, and MEGATEM have the largest footprints given their large transmitter-receiver separations and flight heights. These footprint dimensions are significantly smaller than a typical AEM survey with lateral dimensions in the order of tens of kilometers. Because each footprint is a function of geometry, frequencies, transmitted waveform, receiver time gates, and the 3D conductivity distribution, it is difficult a priori to pick an optimal footprint that can be applied to a an entire survey. Because multiple inversion are always needed to properly assess a final model, the footprint is gra-dually expanded during these multiple inversions until an increase in footprint size does not change the inversion results. The experi-enced user only needs two or three inversions to find the proper footprint size. WB62 Cox et al. Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ SYNTHETIC MODEL STUDY To demonstrate the efficiency of our methodology, we have tested our 3D inversion on a synthetic 3D model of multiple conductors buried with varying dips and depths in a homogenous half-space. The sizes and depths of the body are typical of targets expected in the Abitibi greenstone belt. MEGATEM II step (B-field) response data were synthesized at a 90 Hz base frequency with a 42% duty-cycle half-sine transmitter waveform recording inline and vertical components over five on-time and 15 off-time receiver channels. The transmitter flight height was maintained at a constant 120 m ground clearance with six degrees of pitch. The bird was towed 120 m behind and 50 m below the transmitter. Data were simulated every 20 m along seven alternating flight lines of 100 m separation. As the purpose of this synthetic model was to evaluate resolution, no noise was added to the data. The inline and vertical data were also inverted for layered earth models using AirBeo (Raiche et al., 2007). A 100 ohm-m half-space was used as the starting model for both our 3D inversion and the layered earth inversion. Both inver-sions reduced to misfits less than 2% from the same initial misfit of 23%. Figure 3 presents the same 3D perspective of the actual 3D re-sistivity model and the resistivity model recovered from 3D inver-sion. Figure 4 shows a vertical cross section along the center line for the actual model and the models recovered from 3D inversion and layered earth inversions. It is readily apparent that the 3D inversion recovers each body quite well in their correct location, although the resolution decays at depth. The layered earth inversion failed to ac-curately recover any of the targets, despite obtaining a misfit similar to the 3D inversion. Following industry standard practice, the layered earth inversion results were interpolated to populate a 3D conductivity model for which MEGATEM data were then simu-lated. Despite the layered earth inversion obtaining a misfit less than 2%, 3D modeling of the interpolated 3D conductivity model pro-duced from the layered earth inversions yielded a misfit greater than 19% (Figure 5). The important lesson here is that the quoted misfit for a layered earth inversion is largely irrelevant and has little to no bearing on the quality of any quasi-3D conductivity model recov-ered. To that end, we note that Ellis' (1998) conclusion remains as relevant today as it was over a decade ago: "It is important to note that the degree to which the response of a 1D model can reproduce the observed data is no indication of the one-dimensionality of the local geoelectric model, nor can it be used to infer any relationship between the stitched 1D responses and the true geoelectric earth." CASE STUDY: REID-MAHAFFY TEST SITE, ONTARIO AEM-led exploration for volcanogenic massive sulfide (VMS) deposits in the Abitibi greenstone belt of Ontario and Quebec is often hindered by a conductive near-surface layer. The Reid-Mahaf-fy test site was established by the Ontario Geological Survey as part of Operation Treasure Hunt (OTH). As a condition for awarding surveys in OTH, each contractor was required to acquire AEM data over a common survey area consisting of 16 north-south oriented flight lines of about 5.3 km length, with a line spacing of approxi-mately 200 m. At just 90 line-km, Reid-Mahaffy is a small AEM survey. However, it covers an area of well-known geology that is typically representative of Archean terranes in that it has a mod-erately conductive overburden overlying a resistive basement containing a number of conductive graphites and other bedrock con-ductors. Over the past decade, the site has been flown with multiple AEM systems, including Aerodat, DIGHEM, AEROTEM, VTEM, SPECTREM, GEOTEM, and MEGATEM. The Reid-Mahaffy test site is located in the Abitibi Subprovince, immediately east of the Mattagami River Fault. The area is Figure 3. Perspective of (a) the true 3D resistivity model, and (b) the 3D resistivity model recovered from 3D inversion. Figure 4. Vertical cross section of the MEGATEM synthetic model study. (a) The true 3D resistivity model. Note the bodies are 400 m in strike. (b) Our 3D inversion results. (c) The layered earth inver-sion results using AirBeo. Our 3D inversion and layered earth in-versions commenced from a 100 ohm-m homogeneous half-space starting model, and converged to a misfit less than 2% from an in-itial misfit of 23%. 3D AEM inversion WB63 Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ underlain by Archean (about 2.7 b.a.) mafic to intermediate meta-volcanic rocks in the south, and felsic to intermediate metavolcanic rocks in the north, with a roughly east-west-striking stratigraphy. Narrow horizons of chemical metasedimentary rocks and felsic me-tavolcanic rocks have been mapped, as well as a mafic-to-ultramafic intrusive suite to the southeast. North-northwest-striking Protero-zoic diabase dikes are evident from the aeromagnetic data. Copper and lead-zinc vein/replacement and stratabound VMS mineraliza-tion occur in the immediate vicinity. The world-class Kidd Creek VMS deposit occurs to the southeast of the test site (Ontario Geological Survey, 2000). Drill holes within the test site were aim-ing at conductive targets, and encountered massive sulfides and gra-phite with minor copper and zinc mineralization. AEM data from the Reid-Mahaffy test site have previously been interpreted using a variety of 1D AEM methods (e.g., Witherly et al., 2004; Sattel, 2005; Vallée and Smith, 2007, 2009). In what follows, we compare our 3D inversions for DIGHEM frequency-domain helicopter system, the VTEM-35 time-domain helicopter system, and the MEGATEM II time-domain fixed wing towed bird system. For VTEM and MEGATEM systems, we jointly inverted inline and vertical components for both impulse (dB/dt) and step (B) response data. For running the 3D inversion, our software was parallelized over frequency only. The footprint for the inver-sions varies from 260 m for DIGHEM to 600 m for MEGATEM. Each footprint is found by starting with a small footprint domain and running an inversion. This allows the inversions to run rela-tively fast, and the optimal inversion parameters to be found. As the optimal inversion parameters are found, the footprint size is gra-dually increased. When the results do not change with an increase in footprint size, the footprint is deemed large enough. With experi-ence, this takes only a few attempts. The following results represent the first quantitative comparison of 3D inversions obtained from different AEM systems over the same site. DIGHEM system description and inversion The DIGHEM helicopter system was configured with five oper-ating frequencies: 868, 7025, and 56,374 Hz horizontal coplanar, and 1068 and 4820 Hz vertical coaxial. The transmitter-receiver se-paration was 6.3 m for the highest frequency and 8.0 m for the re-mainder. The survey was flown with a nominal bird height of approximately 32 m. The DIGHEM data were inverted for a 3D conductivity model with 500,000 total cells that were 20 m in the inline direction by 25 m in the crossline direction, and varied from 10 to 20 m in the vertical direction. The model had a total depth of 160 m. The footprint of the DIGHEM system was set at 260 m. The 3D inversion for the DIGHEM data converged to a final misfit of 10% from an initial misfit of 70%. The 3D inversion of the DIGHEM data required just 50 minutes on a Linux workstation using five 2.4 GHz processors and 24 GB RAM. Of the 50 minute run time quoted, 20 minutes were required for Green's tensor precalculations, meaning that the actual inversion itself required just 30 minutes. VTEM system description and inversion The VTEM-35 helicopter system was configured at a 30 Hz base frequency with a 42% duty-cycle polynomial transmitter waveform. The system recorded 34 off-time channels of in-loop inline and ver-tical data. All channels except the very early times were used in the VTEM inversion (channels 1-4), which gave a time range of 0.1 to 7.2 ms after current turnoff. Step (B) responses were calculated from the impulse (dB/dt) responses as per Smith and Annan (2000). The survey was flown with a nominal 35 m ground clear-ance. Inline and vertical components were jointly inverted to 3D conductivity models with 500,000 total cells that were 25 m in the crossline direction, 20 m in the inline direction, and varied from 10 to 70 m in the vertical direction. The model had a total depth of 480 m. The footprint of the VTEM-35 system was set at 400 m. The 3D inversion for the VTEM dB/dt data converged to a final misfit of 53% from an initial misfit of 95%. The VTEM B data con-verged from 95% to 55%. This seemingly high misfit is addressed in the subsequent analysis section. The 3D inversion of the VTEM dB/dt data required three hours on a Linux workstation using eight 2.4 GHz processors and 48 GB RAM. Of the three-hour run time quoted, 1.5 hours were required for Green's tensor precalculations, meaning that the actual inversion itself required just 1.5 hours. Computationally, the only difference between dB/dt and B data oc-curs in the frequency- to time-domain transform. Because the Green's tensors were precomputed for the VTEM dB/dt inversion, they were reused for the VTEM B inversion. This means that our subsequent 3D inversion of the VTEM B data required only 1.5 hours on a Linux workstation using eight 2.4 GHz processors and 48 GB RAM. MEGATEM system description and inversion The MEGATEM II system was configured at a 90 Hz base fre-quency with a 42% duty-cycle half-sine transmitter waveform. The system recorded 20 channels (5 on-time, 15 off-time) of inline and vertical data 125 m behind and 50 m below the transmitter. All channels were used for the inversion, with the last time channel end-ing at 5.5 ms. Step (B) responses were calculated from the impulse (dB/dt) responses as per Smith and Annan (2000). Primary field Figure 5. Synthetic (a) inline and (b) vertical components of ob-served MEGATEM data (solid lines) are shown for different chan-nels. The predicted data from 3D inversion is also shown (þ) and has a misfit less than 2%. The predicted data from 3D modeling of the interpolated layered earth results is shown (Δ) and has a misfit of 19%, even though the layered earth inversion had fit the synthetic data less than 2%. This clearly demonstrates that 1D models are not reliable for quantitative interpretation. WB64 Cox et al. Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ stripping was applied to the impulse and step responses. The survey was flown with a nominal 120 m ground clearance. Inline and ver-tical components were jointly inverted to 3D conductivity models with 500,000 total cells that were 25 m in the crossline direction, 20 m in the inline direction, and varied from 10 to 100 m in the vertical direction. The model had a total depth of 670 m. The foot-print of the MEGATEM II system was set at 600 m. The 3D inver-sion for the MEGATEM dB/dt data converged to a final misfit of 25% from an initial misfit of 70%. The 3D inversion of the MEGA-TEM B data converged to a final misfit of 23% from an initial misfit of 68%. The 3D inversion of the MEGATEM II dB/dt data required three hours on a Linux workstation using eight 2.4 GHz processors and 48 GB RAM. Of the three-hour run time quoted, one hour was re-quired for Green's tensor precalculations, meaning that the actual inversion itself required just two hours. The Green's tensor's were reused for the MEGATEM II B-field inversion, which required only two hours on a Linux workstation using eight 2.4 GHz processors and 48 GB RAM. Analysis Figures 6 and 7 show perspective and volume-rendered images of the 3D resistivity model obtained from 3D inversion of the MEGA-TEM II dB/dt data. Notice that the basement conductors are clearly delineated. In Figure 9, we superimpose the MEGATEM II dB/dt conductivity model from line L50 with by geology as identified by drilling (Ontario Geological Survey, 2000). The conductors identi-fied from drilling in Figure 6 are expanded upon in Table 1.We note there is very good agreement. In particular, the overburden is recov-ered very well. Conductors NM-72-2A-361 and NM-73-2-214 con-tain pyrite concentrated in graphitic black shale sections occurring as nodules, disseminations, and fracture fillings (Ontario Geological Survey, 2000). From Figure 9, these two conductors are interpreted to be the same unit. Conductor NM-73-3-456 is a black chlorite-graphite sediment with pyrite occurring as fine fracture filling and crystals near quartz-calcite veinlets. Figure 8 shows all of L40 of MEGA db/dt x- and z-component data. Note that the peaks in the x-component data align with con-ductors, while the z-component shows a minimum. Figure 10 shows a vertical cross section along line L50 of the conductivity models obtained from the 3D inversion of DIGHEM, VTEM dB/dt, VTEM B, MEGATEM II dB/dt, and MEGATEM II B data. As can be seen, all 3D inversions recover three discrete bedrock conductors beneath the conductive overburden. The four time domain images show a slight southward dip on the northern most conductor (right hand side) and an even smaller northward dip on the central conductor. As expected, the MEGATEM II fixed-wing time-domain system has a deeper depth of investigation than the VTEM helicopter time-domain system, which in turn has a deeper depth of investiga-tion than the DIGHEM helicopter frequency-domain system. As such, DIGHEM is able to recover indications or the tops of some conductors but not resolve them to depth. Figure 6. Perspective image of the 3D resistivity model obtained from 3D inversion of MEGATEM II dB/dt data. Figure 7. Volume rendered image of the 3D resistivity model obtained from 3D inversion of MEGATEM II dB/dt data. Table 1. Conductors identified with drilling. Conductor Depths (m) Description NM-73-2-214 50.0-56.3 15% pyrite, 1% sphalerite, minor chalcopyrite; black shale-andesite tuff (graphite schist - brecciated tuff) NM-73-2A-361 84.3-96.2 20% sulphides, Po 15%, pyrite 3%, chalcopyrite 0.5%, trace sphalerite; mixed black shale (graphitic-chloritic) and fractured andesite tuff NM-73-2A-473 110.5-112.5 50% black shale zone with 10% pyrite-pyrrhotite nodules, disseminations and fracture fillings NM-73-3-456 106.5-116.0 5% pyrite; black shale; black chlorite-graphite sediment 3D AEM inversion WB65 Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ Figure 11 shows the observed and predicted data from the MEGATEM dB/dt channel 11 vertical component and a resistivity cross section from the dB/dt results at 175 m depth. The data com-pare well, and this channel corresponds to a weighted misfit of 25% with our aforementioned weighting scheme. Again, the critical point here is that the quoted misfit for an inversion is largely irrelevant. For example, Figure 12 shows the observed and pre-dicted data from channel 9 of the VTEM dB/dt inline component, Figure 9. Vertical cross section along line L50 of the conductivity model obtained from 3D inversion of the MEGATEM II dB/dt data. Geology, as inferred from drilling, has been superimposed (courtesy of Ontario Geological Survey, 2000). Note that the overburden is mapped, and that distinct subvertical bedrock conductors (graphitic shales containing pyrite) are identified. Figure 10. Vertical cross section along line L50 of theconductivity model obtained from 3D inversion of (a) DIGHEM data, (b) VTEM-35 dB/dt data, (c) VTEM-35 B data, (d) MEGATEM II dB/dt data, and (e) MEGATEM II B data. Figure 8. L40 from the MEGATEM db/dt data. (a) The vertical (z) component data, (b) The inline (x) component, and (c) The vertical cross section of the resistivity. In (a) and (b), the observed data are shown as the solid lines, and the predicted data are shown as the dotted lines. WB66 Cox et al. Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ along with the resistivity cross section at 175 m depth. For clarity, we are only showing the data over the three bedrock conductors shown earlier. Visually, these predicted and observed data compare well, arguably as well as the data shown in Figure 11. However, using the definitions of data weights given above, the weighted mis-fit for this channel is 62%. This is due to the fact that the inline component for the VTEM system is zero for a layered earth and is only sensitive to 3D conductivity variations. This creates many zero crossings, and in our data weighting, small amplitude data have the same importance as larger amplitude data. One may also note that these inline component data cannot be included in any 1D inversion methods as their sensitivity to a layered earth is identically zero in the absence of receiver tilt. Resistivity cross sections and maps of the data show the richness of the geology at the Reid-Mahaffy test site; the complexity of which would be very difficult to honor with qualitative interpreta-tion and impossible with 1D inversion. Figure 11 shows the poten-tial plotting point issues with the towed system geometry. Spatially, the strong responses do not actually correlate with the conductor locations, but are off to the sides. Because different flight directions will produce different coupling with the targets, simple sifting or "lagging" does not correct for this issue. With several electrically Figure 11. Map of data and cross section of the 3D resistivity model from the inversion of the MEGATEM dB/dt data. (a) The observed data from channel 11 of the vertical component of the MEGATEM dB/dt data. (b) The predicted data from the same channel. (c) A horizontal cross section of resistivity at a depth of 175 m recovered from the 3D inversion. Figure 12. Map of data and cross section of the 3D resistivity model from the inversion of the VTEM dB/dt data subset to the area above the three basement conductors. (a) The observed data from channel 7 of the inline component of the VTEM dB/dt data. (b) The predicted data from the same channel. (c) A horizontal cross section of resistivity at a depth of 175 m recovered from the 3D inversion. 3D AEM inversion WB67 Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ close conductors as is the case here, all coupling must be taking into account, which is next to impossible with any methods other than 3D modeling. An oft-cited disadvantage of the towed bird AEM system geo-metry is the asymmetry of the responses with alternating flight line directions, and also the question of where one should place any models or structures. Often, vertical conductors are plotted using the receiver position, and flat-lying conductors are plotted under the transmitter position. This leads to much ambiguity in the actual anomaly position, and arbitrary shifting of data positions to deher-ringbone the data profiles. However, as we have shown, this is not an inherent disadvantage of a towed bird AEM system geometry. This is merely a limitation of the interpretation methods applied to the data. With 3D inversion, all these factors are intrinsically taken into account and the ambiguity in position is removed. For delineating steeply dipping conductors, it is important to in-clude the inline component during 3D inversion as it effectively constrains the vertical component to ensure the correct placement of conductors. For example, for a VTEM survey over a vertical plate, the vertical component exhibits two peaks, and the inline component crosses over between them. Inversion of the vertical component will produce only two conductors, one beneath each peak. Inversion of both the inline and vertical components will pro-duce a single conductor at the correct location. We note that in 1D imaging and inversion methods, the inline component adds little to no value. This reflects an inherent limitation of 1D methods to mod-el the inline component. For example, the inline component of a VTEM system flown over a 1D earth model is identically zero, re-gardless of the structure of that 1D earth model. CONCLUSIONS In this paper, we have presented a paradigm shift for the 3D quan-titative interpretation of AEM surveys. Until recently, quantitative AEM interpretations were based exclusively on conductivity depth imaging or layered earth inversion. We demonstrate how 3D AEM inversion is now a practical consideration for any AEM system, whether frequency- or time-domain, helicopter or fixed-wing. We have presented a model study relevant to the Abitibi greenstone belt, and demonstrated how discrete targets could potentially be re-covered quite accurately to 300 m depth from MEGATEM data. This suggests that potential exists for making new discoveries in the Abitibi greenstone belt by reinterpreting the historic MEGA-TEM data with 3D inversion. Similar conclusions can be drawn elsewhere in the world. With our case study from the Reid-Mahaffy test site, we com-pared results for DIGHEM, VTEM, and MEGATEM systems. As expected, DIGHEM does not detect as deep as the time-domain methods, yet it does best recover the conductive overburden. The VTEM and MEGATEM systems resolve the bedrock conductors quite well, and are still able to image the conductive overburden. The MEGATEM II system, due to its larger dipole moment and transmitter receiver offset, is able to see deeper than the VTEM sys-tem. As expected, the B field response for both systems did produce conductive features at depths slightly greater than the corresponding db/dt response. However, there is much greater difference between systems than between B and dB/dt response within systems. Any of the three systems flown over this area and then inverted in a 3D sense would have resulted in delineation of the bedrock conductors. We note that no conclusion can be made about the accuracy of the 3D models based on their misfits alone. The misfit is but one single number attempting the characterize the quality of fit of an entire inversion, across many channels, components, and geologies, and must be used with caution. Close examination of the resulting 3D models and visual inspection of the predicted and observed data are needed to accurately assess any inversion result. One can often find a misfit metric which will show how "accurate" their inversion is, even if the data fit in reality is poor. In addition, some may use a low final misfit to validate the use of a lower dimensionality in the modeling. In no case is a low misfit any guarantee of an accu-rate model. ACKNOWLEDGMENTS The authors acknowledge TechnoImaging for support of this re-search and for permission to publish.We also acknowledge the con-tributions to 3D AEM inversion by Robert Ellis, Art Raiche, and Fred Sugeng. M. S. Zhdanov acknowledges the support of the University of Utah's Consortium for Electromagnetic Modeling and Inver-sion (CEMI). The Reid-Mahaffy DIGHEM, MEGATEM II, and VTEM-35 data were provided by the Ontario Geological Survey, Fugro Air-borne Surveys Corp., and Geotech Ltd., respectively. We thank Marc A. Vallée and Magdel Combrinck for their assistance regard-ing the MEGATEM and VTEM data, respectively. REFERENCES Ontario Geological Survey, 2000, Airborne magnetic and electromagnetic surveys, Reid-Mahaffy Airborne Geophysical Test Site Survey: Ontario Geological Survey, Miscellaneous Release - Data (MRD)-55. Combrinck, M., L. H. Cox, G. A. Wilson, and M. S. Zhdanov, 2012, 3D VTEM inversion for delineating subvertical shear zones in the West African gold belt: Presented at 22nd ASEG Conference and Exhibition, Brisbane, doi: 10.1071/ASEG2012ab374. Cox, L. H., G. A. Wilson, and M. S. Zhdanov, 2010, 3D inversion of airborne electromagnetic data using a moving footprint: Exploration Geophysics, 41, 250-259, doi: 10.1071/EG10003. Cox, L. H., and M. S. Zhdanov, 2006, Rapid and rigorous 3D inversion of airborne electromagnetic data: 76th Annual International Meeting, SEG, Expanded Abstracts, 795-799, doi: 10.1190/1.2370376. Cox, L. H., and M. S. Zhdanov, 2007, Large scale 3D inversion of HEM data using a moving footprint: 77th Annual International Meeting, SEG, Ex-panded Abstracts, 467-471, doi: 10.1190/1.2792464. Ellis, R. G., 1995, Airborne electromagnetic 3D modelling and inversion: Exploration Geophysics, 26, 138-143, doi: 10.1071/EG995138. Ellis, R. G., 1998, Inversion of airborne electromagnetic data: Exploration Geophysics, 29, 121-127, doi: 10.1071/EG998121. Ellis, R. G., 2002, Electromagnetic inversion using the QMR-FFT fast in-tegral equation method: 72nd Annual International Meeting, SEG, Ex-panded Abstracts, 21-24, doi: 10.1190/1.1817145. Hursán, G., and M. S. Zhdanov, 2002, Contraction integral equation method in three-dimensional electromagnetic modeling: Radio Science, 37, doi: 10.1029/2001RS002513. Leggatt, P. B., P. S. Klinkert, and T. B. Hage, 2000, The Spectrem airborne electromagnetic system - further developments: Geophysics, 65, 1976-1982, doi: 10.1190/1.1444881. Macnae, J., 2007, Developments in broadband airborne electromagnetics in the past decade, in B. Milkereit, ed., Proceedings of Exploration 07: Fifth Decennial International Conference on Mineral Exploration, 387-398. McGillivray, P. R., D. W. Oldenburg, R. G. Ellis, and T. Habashy, 1994, Calculation of sensitivities for the frequency-domain electromagnetic pro-blem: Geophysical Journal International, 116, 1-4, doi: 10.1111/gji.1994 .116.issue-1. Oldenburg, D. W., D. Yang, and E. Haber, 2010, Inversion of multisource TEM data with applications to Mt. Mulligan: 80th Annual International Meeting, SEG, Expanded Abstracts, 3888-3892, doi: 10.1190/1 .3513660. Pare, P., A. V. Gribenko, L. H. Cox, M. Cuma, G. A.Wilson, M. S. Zhdanov, J. Legault, J. Smit, and L. Polome, 2012, 3D inversion of SPECTREM and ZTEM data from the Pebble Cu-Au-Mo porphyry deposit, Alaska: Exploration Geophysics, 43. WB68 Cox et al. Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/ Raiche, A., 1998, Modelling the time-domain response of AEM systems: Exploration Geophysics, 29, 103-106. Raiche, A., F. Sugeng, and G. Wilson, 2007, Practical 3D EM inversion - P223F software suite: presented at ASEG 19th Geophysical Conference and Exhibition, Perth, doi: 10.1071/ASEG2007ab114. Sattel, D., 2005, Inverting airborne electromagnetic (AEM) data using Zoh-dy's method: Geophysics, 70, no. 4, G77-G85, doi: 10.1190/1.1990217. Smith, R. S., and A. P. Annan, 2000, Using an induction coil sensor to indirectly measure the B-field response in the bandwidth of the transient electromagnetic method: Geophysics, 65, 1489-1494, doi: 10.1190/1 .1444837. Smith, R. S., D. Fountain, and M. Allard, 2003, The MEGATEM fixed-wing transient EM system applied to mineral exploration - A discovery case history: First Break, 21, 71-75. Vallée, M. A., and R. S. Smith, 2009, Application of Occam's inversion to airborne time-domain electromagnetics: The Leading Edge, 28, 284-287, doi: 10.1190/1.3104071. Vallée, M. A., and R. S. Smith, 2009, Inversion of airborne time-domain electromagnetic data to a 1D structure using lateral constraints: Near Surface Geophysics, 7, 63-71. Wilson, G. A., L. H. Cox, M. Hope, and M. S. Zhdanov, 2011a, Practical 3D\ AEM inversion - A case study from Golden Ridge, Tanzania: presented at GeoSynthesis. Wilson, G. A., L. H. Cox, and M. S. Zhdanov, 2010, Practical 3D inversion of entire airborne electromagnetic surveys: Preview, 146, 29-33. Wilson, G. A., T. Munday, A. Fitzpatrick, L. H. Cox, and M. S. Zhdanov, 2011b, Practical 3D inversion of AEM data for environmental applica-tions in complex regolith settings: presented at Symposium for Applied Geophysics in Environmental and Engineering Problems, Charleston. Witherly, K., R. Irvine, and M. Godbout, 2004, Reid-Mahaffy Test Site, Ontario, Canada-An example of benchmarking in airborne geophysics: 74th Annual International Meeting, SEG, Expanded Abstracts, 1202- 1205, doi: 10.1190/1.1843294. Xiong, Z., A. Raiche, and F. Sugeng, 1999, A volume-surface integral equation for electromagnetic modeling, in M. L. Oristaglio, and B. R. Spies, eds., Three-dimensional electromagnetics: SEG, 90-100. Zhang, Z., 2003, 3D resistivity mapping of airborne EM data: Geophysics, 68, 1896-1905, doi: 10.1190/1.1635042. Zhdanov, M. S., 2002, Geophysical inverse theory and regularization problems: Elsevier. Zhdanov, M. S., 2009, Geophysical electromagnetic theory and methods: Elsevier. Zhdanov, M. S., and A. Chernyavskiy, 2004, Rapid three-dimensional inversion of multi-transmitter electromagnetic data using the spectral Lanczos decomposition method: Inverse Problems, 20, S233-S256, doi: 10.1088/0266-5611/20/6/S14. Zhdanov, M. S., and E. Tartaras, 2002, Three-dimensional inversion of multi-transmitter electromagnetic data based on the localized quasi-linear approximation: Geophysical Journal International, 148, 506-519, doi: 10 .1046/j.1365-246x.2002.01591.x. 3D AEM inversion WB69 Downloaded 22 Aug 2012 to 155.97.11.184. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/