The deep geological structure of the northwestern United States and southwestern Canada has been extensively studied by seismologists during recent years with the deployment of the National Science Foundation (NSF) EarthScope project’s USArray program, managed by the Incorporated Research Institutions for Seismology (IRIS), and regional seismic networks [1-3]. One of the objectives of these studies was imaging the extent of the subducted Juan de Fuca plate in the mantle east of the Cascades beneath Oregon, and understanding the character of its interaction with the Yellowstone hotspot plume. It was shown in these publications that the subducting Juan de Fuca slab is expressed in the seismological models as a high-velocity anomaly which dips at approximately 46° to the east and extends down to a depth of about 400 km. For example, Obrebski et al. [3], demonstrated that the subducting Juan de Fuca slab is clearly imaged by their P and S-wave tomographic models, which also show a low velocity anomaly associated with the Yellowstone hotspot. Xue and Allen [1] suggested that the absence of the slab below 400 km today is due to the arrival of the Yellowstone plume head about 17 Ma, which destroyed the Juan de Fuca slab at depths greater than the thickness of the continental lithosphere.

An important part of the EarthScope program is the EMScope project, which is the magnetotelluric (MT) component of the USArray program, managed by Oregon State University on behalf of IRIS. By the end of 2011, long-period MT data had been collected at 330 stations located in the following states: Oregon, Washington, Idaho, Wyoming, Montana, California, Nevada, Utah, and Colorado. Preliminary interpretation of the EarthScope MT data collected over Washington, Oregon, Montana, and Idaho were presented by Patro and Egbert [4] and Zhdanov et al. [5]. Zhdanov et al. [11], presented one of the first 3D geoelectrical models of the upper mantle beneath Yellowstone based on the 3D inversion of EarthScope MT data from Montana, Idaho, and Wyoming. These images showed a highly conductive body associated with the tomographically imaged mantle plume-like layer of hot material emerging from the upper mantle toward the Yellowstone volcano. The conductive body identified in those images was westdipping in a similar way to a P-wave low-velocity body [7]. In the recent paper by Kelbert et al. [8], the authors combined spatially sparse EarthScope MT data and higher-resolution MT profiles to produce a geoelectrical model to the depth of 200 km only. This model showed a somewhat surprising result with no conductivity anomaly extending beneath Yellowstone at mantle depths, which contradicted to seismic models.

In the present paper, we have applied our large-scale 3D MT inversion to EarthScope MT data collected over Oregon, Washington, Idaho, Wyoming, Montana, California, Nevada, Utah, and Colorado up to 2011. As a result, we were able to produce a 3D geoelectrical model of the entire northwestern United States, which correlates reasonably well published tomographic models.

Summary of the MT inversion method based on the integral equation approach

In our geoelectrical study of the Earth’s interior beneath the northwestern United States, we employ the MT method, which makes use of variations in natural electromagnetic fields to determine the electrical structure of the Earth. In order to produce a geoelectrical model of the Earth’s crust and upper mantle, one should use a very large-scale 3D inversion covering tens of thousands of square kilometers and extending to hundreds of kilometers in depth. This problem can be solved using massively parallelized 3D MT inversion software. The inversion of MT data represents a classical ill-posed inverse problem; we use regularization theory to obtain a stable solution [9,10]. Considering the extremely large scale of the inversion of the entire Earth Scope MT data collected in nine states of the northwestern United States, we use a moving sensitivity domain (footprint) approach [11]. The interested readers may find more details about our integral equation-based inversion method and of the corresponding large-scale massively parallel inversion algorithm in Zhdanov [12] and Zhdanov et al. [6,11].

Galvanic distortions represent a very well known problem in MT data interpretation. Distortions of the electric field are caused by the accumulation of excess charges formed at high frequency in nearsurface conductivity heterogeneities [13-15]. Over the years, many approaches have been developed to deal with static shift [16,17]. Prior to applying 3D inversion to EarthScope MT data, we established that the effect of the near-surface in homogeneities can be effectively taken care of by 3D inversion, which includes the near-surface inhomogeneous layers. In other words, the 3D inversion models the near-surface in homogeneities well, and thus there is no need for an independent staticshift correction. We have studied carefully the effect of the near-surface in homogeneities on the inversion results for a number of synthetic models, which are not included here due to limited size of the journal paper.

Inversion of Earth Scope MT data

We have analyzed all EarthScope MT data acquired at 330 stations, collected in the northwestern United States by the end of 2011. The station coverage extends from ∼ 123° W to ∼105° W and from ∼38° N to ∼49° N. Original MT data contain four components of the impedance tensor, as well as two components of the magnetic tipper, with error estimations provided for all data points. In our inversion we fit both the amplitude and phase of all four components of the impedance tensors for frequency range between 10s and 10,000s. The considered range of periods corresponds to a depth of investigation down to ∼500 km in a 100 ohm-m host medium. The inversion domain was spanned in the X (geographic E-W), Y (geographic N-S), and Z (vertical downward) directions extending 1650 km, 1350 km, and 500 km, respectively. The sensitivity domain (footprint) size for each station was determined as 450 km based on the rate of sensitivity attenuation of the MT data in a model of a 100 ohm-m half space [11].

Due to the ill-posedness of MT inversion, one can generate a set of equivalent inverse models, which all fit the observed data with the same accuracy [10]. That is why we ran several inversions with different combinations of the data (e.g., full impedance tensor vs. principal impedances), regularization parameters, and discretization of the model, in order to determine the optimal inversion parameters and to examine the robustness of the inversion results. Here, we present just one 3D earth model, which serves as a typical and robust representation of our 3D MT inversion.

The inversions were run on the Ember cluster maintained by the University of Utah’s Center for High Performance Computing (CHPC). We used 48 twelve-core nodes, splitting the work into 96 MPI processes, each of which ran six Open MP threads. The inversions were run until the L2 norm of the residuals between the observed and theoretically predicted MT data, normalized by the L2 norm of the observed data, had decreased to about 10%. The typical inversion runtime was between 12 and 18 hours. For the model presented in this paper, we used cells with a horizontal discretization of 10 km by 10 km, and a vertical discretization starting from 1 km at the surface and logarithmically increasing with depth. The inversion domain contained 2,138,400 cells, and had an initial model of a 100 ohm-m half space.

The conductivity distribution as recovered by our 3D inversion of the EarthScope MT data reflects the regional features of the northwestern United States. Our model correlates well with the seismic P and S-velocity models obtained from inversion of body-wave travel times and surface-wave phase velocities [18]. We present images for different profiles of investigation. Figure 1 shows major geologictectonic feature of northwestern North America overlaid on topography and bathymetry [10]. They encompass the Cascadia subduction zone, including the descending Juan de Fuca and Gorda plates, which are the remnants of the convergent plate margin spanned over coastal western North America. According to modern geological and geophysical data, the Juan de Fuca plate is subducting beneath the North American plate. The estimated depth of the top of the subducting slab is shown with blue contours (labeled in km) in figure 1. The locations of all M>4 earthquakes with depths≥ 35 km since 1970 are shown as blue dots. Volcanoes are shown as orange triangles. The Snake River Plain traces the path of the North American plate over the Yellowstone hotspot, now centered in the Yellowstone National Park (YNP). The Columbia River Flood Basalt Province represents a massive outpouring of basalt from ∼16.6 to ∼15.0 Ma and is shown in pink [19].

Figure 1: Physiographic map of the northwestern United States indicating primary geologic-tectonic features (modified from Obrebski et al., 2011). BM=Blue Mountains; BFZ=Blanco Fracture Zone; BR=Basin and Range; CR=Coast Range; IB=Idaho Batholith; HLP=High Lava Plains; MTJ=Mendocino Triple Junction; YNP=Yellowstone National Park. The locations of EarthScope MT stations used in the inversion are shown by black stars.

Figure 2 presents horizontal sections of geoelectrical model at depths of 100 km, 200 km, and 300 km (panels B, C, and D), respectively. Panel a of figure 2 shows a horizontal section at depth of 100 km overlapped with a map of the major geologic-tectonic features of the Pacific Northwest of the United States, shown in figure 1. The subducting Juan de Fuca slab is clearly imaged in these horizontal sections. It is characterized by the zone of very high resistivity (≥1000 ohm-m) shown by the dark blue vertical strip R1 in panels A, B, C, and D), which corresponds well to the known fact that the subducting oceanic lithosphere is very resistive [20]. We observe several conductive lineaments in the resistivity maps at 100 km depth, similar to those identified by Patro and Egbert [4]. For example, conductive lineament C1 extends under the High Lava Plains to the NW toward the Blue Mountains (panels A and B). It can be associated with the northwestward propagating High Lava Plains volcanic lineament, and it corresponds to the reduced velocity zone (RVZ) beneath the Newberry region, identified by Roth et al. [2], from high-resolution, 3D P-wave tomography. Roth et al. [2], explained the low-velocity anomaly beneath Newberry and the High Lava Plains by the presence of partially melted material due to the release of fluids from the down going slab. This partially melted material may also result in the decreased resistivity (~5 ohm-m), which is observed beneath the western side of the Blue Mountains. This hypothesis is supported by the geochemical characteristics of the Newberry basalts, which, according to Carlson et al. [21], are strongly influenced by subducted slab-derived fluids. Another conductivity anomaly extends beneath the Columbia River basalts in south-central Washington (C2). This anomaly (~5 ohm-m) corresponds well to the reduced velocity zone identified by Roth et al. [2], beneath the thick sections of the Columbia River basalts as related to the geothermal processes in the surrounding areas.

Figure 2: panel A: horizontal section of the geoelectrical model at depth of 100 km, overlapped with a map of the major tectonic features of the northwestern United States. Panel B, C, and D: horizontal sections of the geoelectrical model at depths of 100 km, 200 km, and 300 km, respectively. The dashed contour in B, C, and D denotes the Rocky Mountain front, which serves as the boundary between the tectonically and magmatically active Cordillera to the west and the stable interior of North America to the east.

An important feature of our geoelectrical model is conductive lineament C3 extending beneath Yellowstone-Snake River Plain (YSRP) (Figure 2, panels A and B). This conductive structure corresponds well to the reduced velocity anomaly within the Precambrian lithosphere identified by both Roth et al. [2], and Obrebski et al. [3,18], and is characterized by high temperatures and low density. We observe also conductivity anomaly C4 rising from the mantle at a depth of ∼100–200 km beneath Yellowstone, associated with the Yellowstone conductive plume-like layer identified by Zhdanov et al. [11], using a subset of the EarthScope MT survey of 28 MT stations located over the area surrounding Yellowstone National Park. The highly conductive body is associated with the tomographically imaged mantle plume-like layer emerging from the upper mantle toward the Yellowstone volcano [7,22-27]. Similar to the slow velocity anomaly imaged by James et al. [27], the conductive plume is in fact a mantle layer that extends in a southwest direction into eastern Idaho. This observation opens a possibility for different interpretations of Yellowstone velocity and conductivity anomaly, including tectonic models alternative to “a simple deepmantle plume hypothesis” [3,28].

In horizontal sections of the geoelectrical model (Figure 2), we observe a resistive anomaly (R2) beneath the northern Idaho close to the Idaho Batholith (IB) near the margin of Precambrian North America. Notably, this resistive zone corresponds to the region of increased (~2%) velocities described by Roth et al. [2]. However, contrary to the region of increased velocities observed by Roth et al. [2], which extends to the somewhat unexpected depths of 400 km, the resistive zone extends to depths of up to 200 km, which seems to be more typical for the continental lithosphere. Roth et al. [2] explain the greater than expected depth of the continental lithosphere determined by seismic data by downward smearing of the tomographic image, given the currently limited data coverage in the model for that region.

Figures 3 and 4 show comparisons of the MT inversion result with the DNA09 seismic tomography model of Obrebski et al. [18]. The MT inversion result is presented as the vertical resistivity sections. Vp and Vs seismic velocities are shown as deviations from the mean value in percent. In figure 3, we present a comparison between the vertical sections of the geoelectrical model and the P-wave and S-wave velocity models for the profiles, proceeding along the lines at altitude 40° N, and 42° N, respectively. Figure 4, panel A, shows the similar vertical sections for the profile, proceeding along the line at altitude 44° N. The subducting Juan de Fuca slab is clearly imaged by the geoelectrical model as a high-resistivity anomaly (R1) dipping east under 45°. It is manifested in the P and S-wave tomographic models of Obrebski et al. [13,18] as the high-velocity anomaly which dips east under a similar ∼46°. We also can see the conductive zone of partially melted material directly above the subducting slab, which can be explained by the release of fluids from the down-going slab. Vertical resistivity sections of the geoelectrical model of the northwestern United States show large zones of moderate-to-high conductivity (~5-10 ohm-m) below 100-200 km in the upper mantle, which represent the electrical properties of the conductive electrical asthenosphere (C5).

Figure 3: Comparisons of the MT inversion result with the DNA09 seismic tomography model of Obrebski et. al. (2011). The MT inversion result is presented as the vertical resistivity sections. Vp and Vs seismic velocities are shown as deviations from the mean value in percent. We present the vertical sections of the geoelectrical model and the P-wave and S-wave velocity models for the profiles, proceeding along the lines at latitude 40° N (panel A) and 42° N (panel B).

Figure 4: Comparisons of the MT inversion result with the DNA09 seismic tomography model of Obrebski et. al. (2011) for the profile proceeding along the line at latitude 44° N (panel A), and for profile AA’ along the Snake River Plain (panel B).

Another remarkable geoelectrical feature shown in the vertical resistivity sections (Figure 4) is an extensive area (C6) of low resistivity (~1-10 ohm-m) in the upper mantle, and in some parts, in the low crust, which extends beneath the northwest Basin and Range (BR), High Lava Plains (HLP), Snake River Plain (SRP), and Blue Mountains (BM). Note that, a similar result was observed by Patro and Egbert [4] using 3D inversion of the EarthScope MT data acquired in 2006 and 2007. Interestingly, as early as in 1977, Stanley et al. [29] conducted MT soundings along a profile extending from the Raft River geothermal area in southern Idaho to Yellowstone National Park in Wyoming. The 1D interpretation of these MT sounding curves revealed a highly conductive crustal anomaly with the depth of the conductive zone about 25 km and the resistivity less than 10 0hm-m and at some sites less than 1 0hm-m.

Figure 4, panel B, presents the vertical sections of the geoelectrical model and P-wave and S-wave velocity models for the profile AA’, proceeding along the Snake River Plain (see Figure 1) over the traces of the path of the North American plate over the Yellowstone hotspot. The Yellowstone hot conductive layer is clearly imaged in the geoelectrical model as a low-resistivity (~1-5 ohm-m) anomaly (C4). It is also manifested in the P and S-wave tomographic models of Obrebski et al. [13,18] as the low velocity zone. Note that, in the recent paper by Kelbert et al. [8], the authors combine spatially sparse EarthScope MT data and higher-resolution MT profiles to produce a geoelectrical model to the depth of 200 km only, which shows no conductivity anomaly extending beneath Yellowstone at that depths. As correctly stated by the authors, “This result stands in contrast to the seismic images, which show substantial slow anomalies in the mantle immediately beneath Yellowstone.” This result also contradicts our geoelectrical model which shows the conductive plume-like layer C4 in the mantle (Figure 4, panel B), similar to the models of Zhdanov et al. [11] and the many tomographic models cited in this paper. This drastically different result, published by Kelbert et al. [8], illustrates ambiguity of the MT inversion and significant difficulties with the inversion of long-period MT data, which is similar to the known ambiguity in potential field (gravity and magnetic) inversion. In potential field inversion, a standard approach to avoid unrealistic shallow anomalies is to use appropriate model weights during the inversion which account for attenuation of the sensitivity of the data with the depth [10]. We use a similar technique in our 3D MT inversion method, which reduces the possible ambiguity of the inversion and produces the geoelectrical anomalies located at a depth consistent with seismic images.

The conductivity of the layer C4 is higher in the relatively shallow part (~150 km) and decreases with the depth, diffusing below 300 km. This observation correlates well with the seismic models of Obrebski et al. [13,18], which are characterized by the lower velocity anomaly at the shallow, elongated part of the layer extending to depths of about ~300 km. Both of these observations are consistent with the presence of partial melt, which decreases the seismic velocity and increases the electric conductivity. As was discussed in Zhdanov et al. [11], the low resistivity of the conductive structure (on an order of 1 to 5 ohm-m) is comparable to the resistivity of silicate melts determined from laboratory experiments [30], and it is explained by a combination of high temperature partial melt of basalt and olivine, and the presence of super-critical water present in magmatic processes [31-33]. The low velocities and high conductivity layer initially dip slightly toward the northwest in the upper mantle, consistent with the northwest dip of the mantle plume-like layer seen in published images beneath Yellowstone [5,11].

We have inverted the EarthScope MT data acquired to the end of 2011 over the northwestern United States. Similar to published seismic tomography models, our inverse geoelectrical model of the Earth’s interior beneath the northwestern United States shows a resistive structure associated with the Juan de Fuca slab subducting beneath the Precambrian northwestern United States, and the conductive anomaly characterizing the partially melted material above the subducting slab. The geoelectrical model also contains several prominent conductive features, such as conductive lineaments beneath the High Lava Plains and the Snake River Plain, the conductivity anomaly extending beneath the Columbia River basalts, the conductive mantle layer of the Yellowstone hotspot, and extensive areas of low resistivity in the upper mantle and in the low crust. These results generally correlate well with the P-wave and S-wave velocity models obtained from seismic tomography.

In conclusion, we acknowledge that our geoelectrical model of the northwestern United States represents just one of the first models obtained from 3D inversion of Earth Scope MT data. Such largescale MT inversions are extremely complex, and require significant additional efforts to fully model, understand, and interpret the Earth Scope MT data. At the same time, we believe that the geoelectrical models obtained from our 3D MT inversion provide important complementary information to the published seismic models, and will help to better understand the complex tectonic processes responsible for the formation of the unique geological features of the North America subduction zone. The focus of future research should be on the integrated interpretation of seismic, electromagnetic, gravity, magnetic, and geothermal data, which would reduce ambiguity of geophysical inversion.

We acknowledge the support of the University of Utah’s Consortium for Electromagnetic Modeling and Inversion (CEMI). The MT data were acquired by the Incorporated Research Institutions for Seismology (IRIS) as part of the USArray program of the National Science Foundation (NSF) EarthScope project. Data used in this study were made available through EarthScope (www.earthscope.org; EAR-0323309). We are also thankful to Dr. Richard Allen for providing access to his seismic tomographic models. We acknowledge an allocation of computer time provided by the University of Utah’s Center for High Performance Computing (CHPC). The authors also acknowledge TechnoImaging for support of this research.