Three-dimensional magnetotelluric modeling in the spherical and Cartesian coordinate systems: A comparative study

1.   Introduction

Magnetotelluric forward modeling and inversion modeling are traditionally performed in Cartesian coordinates under the implicit assumption that the modeling domain is sufficiently small such that all numerical considerations caused by the Earth’s curvature may be safely ignored (Kaufman and Keller, 1981; Berdichevsky and Dmitriev, 2008; Chave and Jones, 2012; Liu HY et al., 2022). Indeed, magnetotelluric studies to date have been concerned with sufficiently small spatial scales such that this assumption is easily justified. From radio and audio magnetotelluric applications (Yang B et al., 2019) concerned with regions of a few meters to a few kilometers (Han Q et al., 2021), to wideband and long-period magnetotelluric data intended to probe areas up to a few thousand kilometers in diameter (Xu S et al., 2019), these applications were clearly well within the range of applicability of the Cartesian approximation. The forward methods generally used in the Cartesian magnetotelluric modeling mainly include the integral equation method (Wannamaker, 1991), the finite difference (FD) method (Mackie et al., 1994; Kelbert et al., 2014), the finite volume method (Haber and Ascher, 2001; Jahandari and Farquharson, 2014), and the finite element method (Ren ZY et al., 2013; Cai HZ et al., 2021).

This situation changes with the onset of continental-level magnetotelluric programs, such as EarthScope (Meltzer, 2003), AuScope (Woodcock et al., 2010), and SinoProbe (Dong SW et al., 2013), which are designed to probe the electrical conductivity in the Earth’s interior on a continental scale. In conjunction with a dramatic enhancement in computational capabilities, these programs present an opportunity to apply magnetotelluric techniques to much larger spatial scales than ever before, perhaps reaching well into the Earth’s transition zone with arrays that extend across continents, from shore to shore. These new applications of magnetotellurics require that we revisit the traditional Cartesian approximation and explore the numerical effects of the many implicit approximations that such an approach presents. The major goal of this research was to study the errors that the Cartesian approximation could introduce, relative to global regional magnetotellurics modeling performed in spherical coordinates, and to outline the limits of applicability of the Cartesian approximation in magnetotellurics.

The study of the influence of the Earth’s curvature on electromagnetic responses can be traced back to 1966. Srivastava (1966) derived an expression of impedance for one-dimensional (1D) layered spherical and flat Earths, finding that the Earth’s curvature affects the impedance only for periods greater than one day. In fact, for the 1D electrical conductivity model, a mathematical conversion relationship exists between the spherical and Cartesian coordinates. Weidelt (1972) was the first to propose this conversion equation, later named the “Weidelt transformation.” This method is based on the equivalence of the electromagnetic wave attenuation and the skin depth; these two components are preserved by changing the conductivity value and depth of the model. Conversion from spherical to Cartesian coordinates leads to an increase in conductivity and a decrease in depth. Therefore, the Weidelt transformation can be understood as a type of projection in the vertical direction. Berdichevsky and Zhdanov (1984) and Schmucker (1987) further studied this issue and concluded that the Earth’s curvature is negligible within a depth of 1,400 km. However, these studies focus on 1D models and differ, in terms of spectrum analysis, slightly from the current magnetotelluric method.

In the past few years, the effect of the Earth’s curvature on electromagnetic modeling has regained attention because of the need for continuous large-scale three-dimensional (3D) electrical conductivity models. Numerical simulations are applied to calculate the response of 3D electrical conductivity models in spherical coordinates. Grayver et al. (2019) calculated the impedance tensor in a spherical Earth with a 3D conductivity distribution based on the high-order finite element method and established an equivalent source that results in a valid global tensor. Luo W et al. (2019) studied 3D magnetotelluric forward modeling with a staggered grid in spherical coordinates referring to geomagnetic sounding. Simultaneously, Han Q et al. (2020) used the same grid setting but took the electric field as the primary field to operate the spherical magnetotelluric simulations. All these studies compared the impedance or resistivity difference between traditional Cartesian and spherical coordinate results. The effects of the Earth’s curvature on a 3D electrical conductivity model were found to be much stronger than those on a 1D model and to depend heavily on the geographic projection. However, a number of problems remain to be investigated. For example, large-scale magnetotelluric modeling in spherical coordinates is undoubtedly more accurate, but this does not mean that Cartesian coordinates are no longer applicable. Accordingly, a quantitative standard needs to be discussed and proposed.

In this study, on the basis of a realistic 3D electrical conductivity model, we further explore the errors in magnetotelluric modeling between spherical and Cartesian coordinate models. The three overarching goals of this research are as follows: (1) to present a flow chart of the model conversion from spherical to Cartesian coordinates; (2) to incorporate the usage of meridian convergence in impedance rotation; and (3) to quantitatively assess whether a study area, such as the western or contiguous United States, can still be interpreted using Cartesian coordinates.

2.   Methodology

Here, we first give a brief review of magnetotelluric forward theory. After eliminating the magnetic fields, the 3D quasi-static Maxwell equations for magnetotellurics are written as a second-order elliptic system of partial differential equations in terms of the electric fields alone (using the time dependence of ${\mathrm{e}}^{\mathrm{i}\omega t}$):

where $\boldsymbol{E}$ is the electric field, $\omega$ is the angular frequency, $\mu$ is the magnetic permeability, and $\sigma$ is the electrical conductivity. Additionally, the magnetic fields H corresponding to the electric fields according to Faraday’s law can be expressed as

The magnetotelluric impedance tensor Z maps the horizontal magnetic field ${\boldsymbol{H}}_{{h}}$ onto the horizontal electric field ${\boldsymbol{E}}_{{h}}$ at the Earth’s surface:

The tensor $\boldsymbol{Z}$ is the fundamental measured parameter of magnetotellurics and is the point at which experiment and theory are connected (Tong XZ et al., 2018). In practical applications, individual elements of $\boldsymbol{Z}$ have to be addressed in a specified coordinate system, and the coordinate representation of $\boldsymbol{Z}$ depends on the selection of this system. In some cases, the response tensor Z needs to be rotated to a specific direction via a clockwise angle $\alpha$; the rotation tensor ${\boldsymbol{R}}$ is

Therefore, the rotated response tensor ${\boldsymbol{Z}}'$ is

The impedance tensor can then be used to compute the apparent resistivity $\rho$ and phase $\varPhi$:

To solve Equation (1), an FD method discretized on a staggered grid in spherical coordinates is used. Following the framework of Kelbert et al. (2014) and Zhang H et al. (2019), we developed a regional magnetotelluric modeling code. A sketch of the grid subdivision is shown in Figure 1a. The global model, including both the conductive Earth and the resistive air, is divided into a large number of curved rectangular cells, excepting the areas that are very close to the poles and the Earth’s core. Our computational domain is denoted with red lines in Figure 1a. The electric fields are defined on the edges of every cell bounded by latitude, longitude, and radius, in terms of which Equation (1) is formulated; additionally, the magnetic fields are defined naturally on the cell faces, and the electrical conductivity is defined in each cell center. Specifically, in spherical coordinates, $\boldsymbol{E}$ is decomposed into three components ${E}_{\theta }$, ${E}_{\phi }$, and ${E}_{r}$. Component ${E}_{\theta }$ begins at a low latitude, spans to a high latitude, and points toward the geographic North Pole; component ${E}_{\phi }$ points from west to east; and component ${E}_{r}$ is directed downward toward the Earth’s core (Figure 1b). In Cartesian coordinates, we use the traditional symbols ${E}_{x}$, ${E}_{y}$, and ${E}_{z}$ that point to grid north, east, and down, respectively (Figure 1c). Han Q et al. (2020) derived an FD simulation in spherical coordinates, and they developed 1D layered models and 3D continent–ocean models to prove the correctness and stability of the regional spherical magnetotelluric forward program.

Figure  1.  (a) Sketch of the grid subdivision in spherical coordinates. (b) Curved rectangular cell in spherical coordinates. (c) Rectangular cell in Cartesian coordinates (modified from Zhang H et al., 2019, and Han Q et al., 2020).

DownLoad: Full-Size Img PowerPoint

3.   Comparative Study

The workflow of magnetotelluric inversion begins with mapping the area of interest and the site locations to some sort of Cartesian framework and ends with the co-location of the Cartesian grid with geographic coordinates for geologic interpretation. In the magnetotelluric practitioner’s toolbox, no standard tools are available to perform these tasks, and a variety of homegrown utilities are used. Formally, these tasks may be described as a model conversion that defines the locations of data points in the artificial Cartesian coordinates. In this study, we explore the effects of model conversion on the accuracy of the forward-modeled responses, as compared with direct modeling in spherical coordinates.

3.1   Model Conversion

The model conversion process is represented in Figure 2a. According to the geographical location information, it is convenient to set up an electrical conductivity model in spherical coordinates that is defined by the latitude, longitude, and radius. This spherical model is then divided into two parts: the spherical electrical conductivity, σ_s, and a spherical grid, grid_s. A map projection method is chosen to project grid_s to Cartesian coordinates, that is, the geography is converted into length in this step. However, because of the use of a regular grid in the FD method, further equalization is necessary to obtain a Cartesian grid, grid_c, as shown in Figure 2b. During the conversion from grid_s to grid_c, the number of grid cells remains the same such that, in terms of mathematics, σ_s can be copied directly to grid_c. Assigning electrical conductivity in this way works only for specific map projections, such as cylindrical equal distance. In most map projection methods, directly copying the electrical conductivity will lead to a serious offset between the magnetotelluric sites and the underground structure. Therefore, we project the center points of the Cartesian cells back to spherical coordinates by using the same map projection method, obtaining a new spherical grid, grid_s2. Combining the known grid_s and σ_s, we can obtain the electrical conductivity σ_c of the Cartesian model by using interpolation methods; the nearest neighbor is used in this paper. The values of σ_c and grid_c then form the Cartesian model.

Figure  2.  (a) Flowchart of the model conversion from spherical to Cartesian coordinates. (b) Sketch maps of a Cartesian grid generated from five representative map projections. The red letters are abbreviations of the map projection methods, and the blue lines or stars refer to nondistortion locations. The characteristics of these five map projections are given in Table 1.

DownLoad: Full-Size Img PowerPoint

Five common map projections (whose characteristics and definitions are described in Table 1 and Appendix A) are used to convert the electrical conductivity model from spherical to Cartesian coordinates. According to the different properties of these map projections, we need to choose different parameters presenting the nondistortion locations, marked by the blue lines or stars in Figure 2b. To maintain the balance of the model, the middle latitude is considered the standard parallel for the cylindrical equal distance (eqdcylin) and cylindrical equal area (eqacylin) projections, and the center point of the model is used for the Universal Transverse Mercator (UTM) and azimuthal equal area (eqaazim) projections as the nondistortion point. The Lambert conformal conic (lambertstd) projection requires two selected standard parallels. Because there are many choices, we use the middle latitude minus and plus one fixed number, keeping these two parallels at one quarter and three quarters of the model. The same map projection is then used to compute the magnetotelluric site locations in the Cartesian coordinate system.

Table  1.  List of the five map projection methods (Snyder, 1987).

Abbreviation	Type	Properties
eqdcylin	Cylindrical	Equal distance, scale is true along all meridians and the standard parallel, scale is constant along any other parallels, distortion of both shape and area increase with distance from the standard parallel.
eqacylin	Cylindrical	Equal area, shape distortion increases with distance from the standard parallel.
UTM	Cylindrical	Conformal, Universal Transverse Mercator system divides the Earth into zones (each 8° × 6° in extent) that use formulas for a transverse version of the Mercator projection, with projection and ellipsoid parameters designed to limit distortion.
lambertstd	Conic	Conformal, scale is true along the two selected standard parallels, distortion is constant along any other parallel, conformal everywhere but the poles, not equal area or equal distance.
eqaazim	Azimuthal	Equal area, only the center point is free of distortion, scale is true only at the center point, increasing tangentially and decreasing radially with distance from the center point.

 | Show Table

DownLoad: CSV

The western United States model (hereafter referred to as WUS) shown in Figure 3 is the original spherical model used for the model conversion and is the result of a 3D magnetotelluric data inversion providing a regional-scale view of the electrical conductivity from the middle crust to nearly the mantle transition zone, covering an area from northwest Washington to northwest Colorado (Meqbel et al., 2014). This WUS model arches across approximately 14 latitudes and 25 longitudes, including part of the Pacific Ocean. The entire area is discretized into 124 × 156 × 43 cells, with the conductivity defined in the center of each cell. The cell thicknesses are 500 m in the upper part of the model, increasing logarithmically downward to a depth of 1,468 km. The full magnetotelluric impedance tensor is sampled at 325 sites distributed on a quasi-regular two-dimensional array with site spacings of approximately 70 km in both horizontal directions.

Figure  3.  Depth section at 1 km of the western United States (WUS) conductivity model. The black dashed lines represent the coastline, and the solid black lines represent the boundaries of the states.

DownLoad: Full-Size Img PowerPoint

Figure 4 shows the model conversion process and the results of the five map projection methods. For each projection, three panels are shown to evaluate its effect from different perspectives. The first panel shows an interpolation diagram, which is a key step showing how to obtain σ_c for the Cartesian model. The red dots represent the locations of the spherical conductivity σ_s, and the blue dots represent the cell centers of grid_s2; that is, the blue dots show the locations where interpolation is needed. The best situation occurs when the blue dots and red dots can completely overlap, such as with the eqdcylin projection. If the blue dots expand outward beyond the red dots, this part will not be interpolated and the corresponding regions of the Cartesian model will be null. If the blue points shrink, they will have repeated values. In other words, the original spherical electrical conductivity indicated by the red dots cannot be accurately passed to the Cartesian model when the blue dots expand or shrink because of the distortion of the map projection. For example, the blue dots from the UTM, lambertstd, and eqaazim map projections present expansion at high latitudes and shrinkage at low latitudes.

Figure  4.  Model conversion results for the five map projections, with the projection method corresponding to the numbers as follows: (a) eqdcylin; (b) eqacylin; (c) UTM; (d) lambertstd; and (e) eqaazim.

DownLoad: Full-Size Img PowerPoint

The second panel shows the transformed Cartesian electrical conductivity model, with the coastline and state boundaries projected using the same method. Visually, these five Cartesian models can be approximately divided into two types according to their deformation: the UTM, lambertstd, and eqaazim projections maintain the curvature of the Earth, whereas the eqdcylin and eqacylin projections seem suitable for Cartesian coordinates and show a more balanced shape. The third panel shows the model difference, which is the difference in the conductivity between the WUS model and its related Cartesian model, taking ${\sigma }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}= \left|{\mathrm{log}}_{10}(\sigma \_s)-{\mathrm{log}}_{10}(\sigma \_c)\right|$. The distribution of the model differences is closely related to the characteristics of the adopted map projection and has an intuitive correlation with the interpolation diagram. We can conclude that basically no distortion occurs in the obtained Cartesian model near the projection standard parallels or centers. However, with increasing distance from the projection standard parallels or centers, the distortion of the Cartesian model becomes increasingly serious.

3.2   Differences Between the Apparent Resistivities

In magnetotelluric sounding, the apparent resistivity is one of the intuitive parameters reflecting the underground electrical conductivity distribution and is often used as an input parameter for inversion. In this section, we further compare the difference between the apparent resistivity calculated in spherical coordinates and that calculated in Cartesian coordinates. To make a quantitative comparison, Equation (8) gives an expression of the apparent resistivity difference ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$:

where ${\rho }_{c}$ is the apparent resistivity calculated in Cartesian coordinates, consisting of two components ${\rho }_{xy}$ and ${\rho }_{yx}$, and ${\rho }_{s}$ is the apparent resistivity calculated in spherical coordinates; its two components that correspond to ${\rho }_{xy}$ and ${\rho }_{yx}$ are ${\rho }_{\theta \phi }$ and ${\rho }_{\phi \theta }$, respectively. For convenience in the following description, we define

The distribution of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ is represented in Figure 5, where ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}-s}$ in each subplot title indicates the average difference over all sites for one period. The results show that the ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ values are the smallest between the Cartesian model projected by eqdcylin and the original spherical model. As introduced in Table 1, the cylindrical equal distance method projects the global model into a cylinder, with the distances along the middle latitude being conserved. Therefore, large ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ values are primarily distributed at relatively high and low latitude regions, and in the area along the middle latitude, namely the middle belt of the model, the difference is nearly zero.

Figure  5.  Apparent resistivity difference between the original spherical model and the transformed Cartesian models in Figure 4 for a period of 100 s. The black dots indicate the locations of the magnetotelluric sites.

DownLoad: Full-Size Img PowerPoint

4.   Meridian Convergence

Although less obvious than the map projection, accurate regional data interpretation requires that all data be oriented along the north of the model grid. For field data that are oriented to geographic north, this requirement holds in spherical coordinates. However, once the model is converted to Cartesian coordinates according to some map projections, this data orientation assumption no longer strictly holds. At all longitudes except for the longitude of the grid center, the field data are now oriented at an angle to the grid north (or Cartesian x-axis): in the Northern Hemisphere, they point inward, whereas in the Southern Hemisphere, they point outward relative to the Cartesian grid north. This phenomenon is best known as meridian convergence and is well recognized and corrected for in geodesy and cartography (Reilly and Bibby, 1975; Soler and Fury, 2000). For continental-scale magnetotelluric applications, meridian convergence is yet another effect contributing to the inaccuracy of the Cartesian approximation and needs to be accounted for and corrected within the framework of this comparison. That is, when we compare the differences in the forward data, we need to ensure that they are in the same direction. The data orientation correction angles can be analytically computed for most map projections. Here, we derive these expressions for the five map projections outlined in Table 1. We then analyze the effect of the meridian correction on the realistic WUS conductivity model.

A schematic diagram of the meridian convergence in the Northern Hemisphere is shown in Figure 6. The computed impedance tensor from the Cartesian model at location P points to grid north, indicated by the straight red arrow, whereas the direction of the spherical impedance is consistent with the tangent red arrow T. If the meridian convergence angle $\gamma$ is known, the Cartesian impedance can be rotated to the T direction. According to the definition, the meridian convergence angle $\gamma$ is zero at the central meridian, is positive when P is east of the central meridian, and is negative when P is west of the central meridian.

Figure  6.  Schematic diagram of meridian convergence. Point P is a given point in the Northern Hemisphere, the central meridian appears as a straight line (the grid north axis), line T is the tangent line to the meridian at P, and $\boldsymbol{\gamma }$ is the meridian convergence angle, the angle measured clockwise from the tangent line to grid north, with a range of [−90°, 90°].

DownLoad: Full-Size Img PowerPoint

Even though many formulas exist for calculating the meridian convergence angle $\gamma$, it can be derived from a single principle. Assuming the mapping equations for P are

where $\varphi$ and $\lambda$ are the latitude and longitude of P, respectively, $f$ and $g$ represent the projection functions mapping latitude and longitude, respectively, to the x–y-axis. The total differentials of x and y are

For the meridian, we can take $\lambda$ as a constant and the total differentials can be rewritten as

Finally, the meridian convergence angle $\gamma$ defined above is given by

As long as the projection function is determined, according to Equation (15), its corresponding meridian convergence angle $\gamma$ can easily be derived. The eqdcylin and eqacylin projections have a convergence angle of 0 because their mapping equation $f$ is not a function of $\varphi$, such that $\dfrac{\partial x}{\partial \varphi }$ is always 0. For UTM,

where ${\lambda }_{0}$ is the center meridian. For lambertstd,

where $n$ is a constant coefficient that is defined by ${\varphi }_{1}$ and ${\varphi }_{2}$; ${\varphi }_{1}$ and ${\varphi }_{2}$ are the two standard parallels or latitudes of the lambertstd map projection. For eqaazim, the tangent of $\gamma$ is a mixed operation of trigonometric functions and can be summarized as

where ${\varphi }_{1}$ is the standard parallel or latitude used in the eqaazim projection.

Equations (16) to (19) indicate that the meridian convergence angle is zero for regular cylindrical projections and is a simple function of latitude for regular conic and polar azimuthal projections. The definitions of other types of map projections can result in more complicated sets of equations. In this case, the meridian convergence angle needs to be carefully derived from the original principle. Appendix A gives the definitions of the five map projections in this study, showing the forward projection from spherical coordinates to Cartesian coordinates and the details of the derivation process for their corresponding meridian convergence angles.

After the meridian convergence angle $\gamma$ is determined, the rotated impedance can be calculated by using Equation (5), and the meridian convergence angle $\gamma$ at the magnetotelluric site is taken to be the rotation angle $\alpha$. We applied this meridian convergence correction to the impedance calculated from the Cartesian models (i.e., rotating Z from grid north to geographic north to maintain consistency with the direction of impedance calculated from the spherical model) projected by the lambertstd, UTM, and eqaazim projections.

To avoid the contingency of the results, we recalculated the forward results of the spherical and Cartesian models by using the same frequency as the field data (from 7.31 to 18,724 s, 30 frequencies in total). Similarly, we then calculated the apparent resistivity and phase for comparison. Figure 7 shows the variation of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}-s}$ with respect to the period. The red lines in Figure 7 represent the original error, and the black lines represent the error after the meridian convergence correction. The results shown in Figure 7 prove that the meridian convergence correction can further reduce the forward difference between the spherical and Cartesian models.

Figure  7.  Variation in ${{\rho }}_{\rm{diff}-\mathit{s}}$ with respect to the periods between the spherical and Cartesian models transformed by the UTM, lambertstd, and eqaazim map projections. The black lines show the error after the meridian convergence correction, and the red lines indicate the error without the correction.

DownLoad: Full-Size Img PowerPoint

To better understand the influence of the meridian convergence correction on the forward results for a single period and magnetotelluric site, we take the lambertstd map projection as an example. Figure 8 shows the distribution of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ in two selected periods, with the base map showing ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ after the meridian convergence correction. All the arrows on the magnetotelluric sites point to geographic north, and the length of each arrow indicates the change in ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ from before to after the meridian convergence correction. A red arrow indicates that ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ has increased, and a green arrow indicates that ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ has decreased. As shown in Figure 8, most sites are marked with green arrows, and the larger amplitude changes are primarily concentrated at lateral sites. Even though ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}-s}$ shows a decreasing trend, the ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ values of some sites marked by red arrows increased after the meridian convergence correction. This situation occurs when ${\rho }_{s}$ matches ${\rho }_{c}$ well and the impedance rotation increases their difference.

Figure  8.  Distribution of ${{\rho }}_{\rm{diff}}$ (spherical model vs. the lambertstd Cartesian model) after the meridian convergence correction for two periods (T = 102.4 s and T = 11,915.64 s). The length of the arrow drawn at each site indicates the change in ${{\rho }}_{\rm{diff}}$, where green indicates that the Cartesian apparent resistivity is closer to the spherical apparent resistivity after the meridian convergence correction, and red indicates the opposite. All the arrows point to geographic north, and their angle with respect to grid north indicates $\mathit{\gamma }$.

DownLoad: Full-Size Img PowerPoint

Figure 9 presents the apparent resistivities and phases of four magnetotelluric sites. The blue circles and red crosses represent the apparent resistivities of the spherical model and the lambertstd Cartesian model, respectively, and the black triangles represent the apparent resistivity after the meridian convergence correction. If the original Cartesian apparent resistivity was a poor fit with that of the spherical model, an obvious correction occurred after the impedance rotation by the angle $\gamma$ (e.g., WYK17, ORG03, and MTG18). These sites are usually located on the sides of the model far from the middle longitude, that is, they have large values of $\gamma$. The situation for site CAQ01 is relatively complicated: the ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ value after the meridian convergence correction increased close to 100 s. However, from either a qualitative or quantitative view, the corrected apparent resistivity fit better with that of the spherical model.

Figure  9.  Apparent resistivities and phases of four magnetotelluric sites. The blue circles and red crosses represent the apparent resistivities of the spherical and lambertstd Cartesian models, respectively, and the black triangles represent the apparent resistivities and phases after the meridian convergence correction.

DownLoad: Full-Size Img PowerPoint

5.   Quantitative Evaluation

In the previous sections, we explained how to perform a model conversion and compared the differences in the apparent resistivities calculated in different coordinate systems.

In this section, we quantitatively evaluate whether Cartesian coordinates are suitable for the WUS model through a comparison with field data. Magnetotelluric inversion algorithms, such as nonlinear conjugate gradient algorithms, usually set a floor error to control the offset between the forward data and the field data (Egbert and Kelbert, 2012). When the offset is within this floor error, the inversion model is considered acceptable. We applied this setting to our quantitative evaluation, that is, if ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ was smaller than the floor error, we concluded that the influence of the Earth’s curvature could be ignored, otherwise it could not be ignored. The floor error is not a fixed value but is usually within 10% in magnetotelluric inversions; here, we used 5% and a calculation formula of $5\mathrm{\%}\times \sqrt[2]{\left|{Z}_{xy}\right|\times \left|{Z}_{yx}\right|}$, where ${Z}_{xy}$ and ${Z}_{yx}$ are the impedances of the field data.

Figure 10 presents the apparent resistivities and phases of the field data. (The field data are shown with error bars, where blue error bars indicate the real measurement error and orange error bars indicate the 5% floor error.) Cartesian and spherical forward data of three representative sites were selected from the 325 magnetotelluric sites. According to the previous comparison results in Section 3, the Cartesian forward data compared here were calculated from the Cartesian model transformed from the eqdcylin projection.

Figure  10.  Apparent resistivities and phases of the western United States model. The blue circles, red crosses, and green stars represent the magnetotelluric field data, spherical forward data, and Cartesian forward data, respectively, the blue error bars represent the real measurement error, and the orange error bars represent the 5% floor error.

DownLoad: Full-Size Img PowerPoint

Figure 10 shows that sites MTB14 and MTD17 have some common characteristics: both the spherical and Cartesian forward data are a poor fit with the field data, and their apparent resistivity curves are obviously separate when compared with site IDI11. However, for MTB14, ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ is larger than the real error but smaller than the 5% floor error, and for MTD17, it is larger than the 5% floor error within a certain period. According to our statistics of the 325 magnetotelluric sites, only 6% of sites are similar to MTD17. The spherical and Cartesian apparent resistivity curves for more than 80% of the sites are coincident with IDI11, fit the field data well, and show almost no difference. In this case, we concluded that the WUS model could still be simulated by using Cartesian coordinates, given an appropriate model conversion method.

6.   Continental-Scale Three-Dimensional Magnetotelluric Comparison

A study area similar to the WUS can still be simulated, inverted, or interpreted by using the Cartesian coordinate system after adopting appropriate map projection methods. If the scope of the study area is further expanded, such as to the contiguous United States, whether Cartesian coordinates are still applicable needs to be reevaluated. Figure 11 shows two depth sections from a 3D electrical conductivity model of the contiguous United States, having a resolution of 0.25° in latitude and 0.5° in longitude (Kelbert et al., 2019). Incorporating the findings from the previous sections, we transformed this continental United States model to its corresponding Cartesian model by using a cylindrical equal distance map projection with a middle latitude of 37°N as the standard parallel, resulting in a Cartesian model with a resolution of 28 × 44 km.

Figure  11.  Depth sections from the three-dimensional electrical conductivity model of the contiguous United States at 10 km (left) and 50 km (right).

DownLoad: Full-Size Img PowerPoint

The apparent resistivities of the total 933 magnetotelluric sites were calculated in the different coordinate systems. Figure 12 demonstrates the spatial distribution of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}-p}$, where ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}-p}$ is the average value of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ over 30 periods for one magnetotelluric site. We concluded that the distribution of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}-p}$ conforms to the rules of a cylindrical equal distance map projection. The value of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}-p}$ is smallest along the middle latitude area and is relatively larger along the east coast and in the northern regions of the contiguous United States.

Figure  12.  Distribution of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}-p}$ for the three-dimensional contiguous United States electrical conductivity model. The magnetotelluric sites are represented by triangles, with the color fill in each triangle indicating the value of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}-p}$.

DownLoad: Full-Size Img PowerPoint

Figure 13 shows the apparent resistivities and phases of three selected magnetotelluric sites (WAC08, MNC39, and WAD10) whose ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}-p}$ values are larger than 0.07, where the separation of the spherical and Cartesian apparent resistivity curves is very obvious.

Figure  13.  Apparent resistivities and phases calculated from the contiguous United States model. The blue circles, red crosses, and green stars represent the magnetotelluric field data, spherical forward data, and Cartesian forward data, respectively, the blue error bars represent the real measurement error, and the orange error bars represent the 5% floor error.

DownLoad: Full-Size Img PowerPoint

Following the quantitative rules of the WUS model, we made corresponding quantitative evaluations for the contiguous United States model. The magnetotelluric sites with ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ values greater than the 5% floor error are marked as black triangles in Figure 14. These magnetotelluric sites are seriously affected by the Earth’s curvature and account for 12% of the total sites. Accordingly, this contiguous United States model may not be suitable for Cartesian coordinates for two main reasons. One is that the distribution of the magnetotelluric sites in Figure 14 is no longer scattered but shows continuity and regionality. This may mislead magnetotelluric interpretations in such areas, especially interpretations that need to use the inverted resistivity to infer the water content, temperature, or both. The second reason is that the contiguous United States model in this study is incomplete; the electrical conductivity in the southern part of the model has extremely low resolution, and the magnetotelluric sites do not cover the entire model area. We predict that the number or percentage of magnetotelluric sites in Figure 14 will increase once the contiguous United States model is complete.

Figure  14.  Distribution of magnetotelluric sites with ${{\rho }}_{\rm{diff}}$ values greater than the 5% floor error. Each site is shown as a black triangle.

DownLoad: Full-Size Img PowerPoint

7.   Discussion

A series of comparative experiments were performed to study the difference between the apparent resistivities calculated in spherical and Cartesian coordinates. Both calculations used numerical simulation methods. Even though the control parameters in the numerical simulations were the same, calculation errors were still inevitable. This led to the following problem: How can it be determined whether the compared difference is caused by the different coordinate systems or by the calculation itself? Consequently, we rarely used a single data point for a comparative analysis but rather compared the average values of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ and examined the overall differences. Strong evidence exists that the distribution of ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ in Figure 5 matches the characteristics of the corresponding map projections, indirectly proving that the influence of the Earth’s curvature is primary.

The comparative experiments in Section 3 showed that the apparent resistivity calculated from the Cartesian model obtained by the eqdcylin projection is the closest to the calculation result of the spherical model. The reason this cylindrical equal distance method performs best may be related to the 3D magnetotelluric simulation method. Both spherical and Cartesian coordinates use the FD simulation method. This method calculates the edge length as a basic component to obtain other geometric elements. The smallest apparent resistivity difference ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ may be generated by the projection method that can maintain the distance or the length without distortion to the greatest extent. The second smallest ${\rho }_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}$ arises from the eqacylin map projection. Note that both eqdcylin and eqacylin are cylindrical projection methods. A cylindrical surface can be expanded into a rectangular plane, which naturally adapts to Cartesian coordinates.

We cannot conclude, however, that the cylindrical projection method is superior to other methods. The characteristics of the selected research area WUS represent only one of the important influencing factors. As shown in Figure 3, the WUS model is located in the mid-latitude zone and primarily extends east–west, which guarantees that the upper and lower boundaries are not far from the middle latitude. If the research area is oriented south–north, or close to the Earth’s poles where the arc length changes significantly, the error of cylindrical projections will increase. In this case, other types of map projection methods, such as UTM or eqaazim, may be more appropriate.

The inversion control parameter, namely the 5% floor error, was used as a benchmark in the quantitative evaluations. The floor error is not a fixed value and can be adjusted according to the quality of the field data. If the floor error is set to a large value (e.g., 10%, which is not rigorous), it may affect the conclusion concerning whether Cartesian coordinates can be used in the contiguous United States model. One obvious advantage of spherical magnetotelluric simulations is that they no longer require the map projection step and can effectively avoid errors caused by improper projections. Moreover, in the subsequent interpretation, the spherical electrical conductivity model matches the surface topography and other geological information well.

8.   Conclusions

When implementing a large-scale magnetotelluric simulation, the traditional Cartesian coordinate approximation may not be suitable because of the nonnegligible curvature of the Earth. On the basis of a real 3D conductivity model, we used a newly developed regional spherical magnetotelluric program to explore the difference between spherical and Cartesian simulations. The apparent resistivity calculated from a Cartesian model obtained by a cylindrical equal distance projection is closest to the results calculated directly from the spherical model. For Cartesian models obtained by noncylindrical projections, their impedance requires meridian convergence correction, which can help reduce their difference from the spherical results. Quantitative evaluations showed that Cartesian coordinates can still be used in the WUS model but that the contiguous United States model is no longer suitable for Cartesian coordinates. At present, the calculation time of the regional spherical magnetotelluric program is approximately twice that of Cartesian programs. Improving the calculation efficiency and completing the inversion code will be a topic of future work.

Appendix A: Definitions of five map projection methods and derivation of the meridian convergence angle

Five widely used map projections were introduced in this article. Here, we present their formulas and describe some critical parameters that control the projections. For the cylindrical equal distance projection,

For the cylindrical equal area projection,

For the Lambert conformal conic projection,

of which

For the Lambert azimuthal equal area projection,

where

The Universal Transverse Mercator (UTM) system is not a single map projection. It divides the Earth into 60 zones, with each zone being a 6° band of longitude, and uses a secant transverse Mercator projection in each zone:

The symbols in the equations above are listed here:

$\lambda$ : the longitude of the location to project.

$\varphi$ : the latitude of the location to project.

${\varphi }_{1}$, ${\varphi }_{2}$: standard parallels where the scale of the projection is true.

${\lambda }_{0}$ : the central meridian of the map.

${\varphi }_{0}$ : the central parallel of the map.

$x$ : the horizontal coordinate of the projected location on the map.

$y$ : the vertical coordinate of the projected location on the map.

$R$ : the radius of the sphere.

${k}_{0}$ : the scale factor along the central meridian.

According to Equation (18), we can deduce the meridian convergence angle $\gamma$ when we know the projection relationship. The derivation of x with respect to $\varphi$ is zero for the cylindrical equal distance and area projections, resulting in $\gamma$ being zero. For the Lambert conformal conic projection, n, ${\tau }_{0}$, and $F$ are constant once the standard parallels are decided:

Therefore,

For UTM, by first replacing $\mathrm{tan}\varphi /\mathrm{cos}\left(\lambda -{\lambda }_{0}\right)$ in Equation (A15) with M, we get

where

Substituting Equations (A22) and (A23) into Equations (A21) and (A20), we get the meridian convergence $\gamma$ of UTM after algebraic operations:

The derivation of $\gamma$ for the Lambert azimuthal equal area projection is tricky. Variables x, y and K are functions of $\varphi$. Here, we set some intermediate variable to simplify its expression, defining the denominator of K as

and the factor of Equation (A12) as

Hence, $\gamma$ can be roughly written in the form of

where $\dfrac{\partial K}{\partial \varphi }=\dfrac{\partial K}{\partial m}\cdot \dfrac{\partial m}{\partial \varphi }$ and $\dfrac{\partial g}{\partial \varphi }$ can easily be computed, then substituting them into Equation (A27) to get $\gamma$.