Mantle Plume Reconstruction by Three-Dimensional Electromagnetic Induction

The Earth’s interior consists of multiscale structures that range from micrometer-scale mineral assemblages to 1000-km-scale heterogeneities. Mantle plumes are one such mega-scale structure that connects the core–mantle boundary with Earth’s surface. Reconstructing these structures can provide insights into mantle material and energy convection, as well as Earth’s long-term evolution. However, mantle plume has not yet been convincingly reconstructed by electromagnetic (EM) induction; even they have significantly high electrical conductivity compared with the surrounding mantle. Here, we numerically reconstruct mantle plumes by employing a deep Earth EM induction method—geomagnetic depth sounding (GDS). We build the electrical structure of mantle plumes and conduct inversion tests to investigate how different station coverage areas, station spacings, noise levels, and response period ranges influence the construction. The test results indicate that the reconstruction of a broad 10° diameter plume head near the mantle transition zone (MTZ) requires a station coverage area of at least 10 $^{\circ }\,\,\times10^{\circ }$ and a 2° station spacing; the station spacing can be increased to 5° for a 20 $^{\circ }\,\,\times20^{\circ }$ coverage area. A continuous two-year record with ~5% noise is sufficient to recover the electrical structure of the plume head. Plumes with different types and roots can be distinguished by images near the MTZ, while reconstruction of the narrow tail in the deep lower mantle seems to be difficult due to the limited resolution of GDS. A mantle plume beneath South China is discovered by GDS from the field geomagnetic data. GDS is expected to be used for reconstructing the mantle plume beneath important locations and contributing to the study of Earth’s dynamics.


I. INTRODUCTION
E LECTROMAGNETIC (EM) induction is widely used to elucidate the structure and properties of the Earth's interior. Unexploded ordnance and buried structures near the surface can be identified by EM induction based on their magnetic susceptibility and polarizability [1], [2], [3], [4], [5]. Groundwater, hydrocarbons, and mineral resources can also be detected based on their associated EM induction responses in the frequency and time domains via airborne [6], [7], controlled-source [8], [9], [10], and transient [11], [12], [13] EM methods. Deeper investigations of geological structures, such as magma chambers beneath volcanoes, hotspot pathways, and other lithospheric and upper mantle structures, require much lower frequencies (i.e., 10 −4 -10 2 Hz) for effective EM imaging, with the magnetotelluric (MT) method serving as a representative EM approach at these lower frequencies [14], [15], [16]. Mantle plumes that originate from the core-mantle boundary (CMB) can transport lower mantle material and energy to the surface. These plumes may be closely related to paleocontinent breakup, the formation of large igneous provinces, extensive intraplate volcanism, and even mass extinction events [17], [18]. Therefore, mantle plume reconstruction plays a key role in understanding mantle convection and Earth's evolution. The low viscosity of thermal mantle plumes makes them rise rapidly in the lower mantle, forming a narrow ∼100-km-diameter tail [19]. However, numerical simulations indicate that the diameter of a mantle plume in the lower mantle may be much broader than previously thought [20]. The upwelling mantle plume will react with layer interfaces in the Earth's interior [21], such as the 660-km mineral phase transition interface and the lithosphere-asthenosphere boundary (LAB). When a plume impinges on the mantle transition zone (MTZ), the exothermic reaction caused by the phase transformation at the 660-km interface can produce a tremendous volume of hot plume material that accumulates beneath the MTZ [21], with the extent of this plume head potentially exceeding 1000 km in diameter [19].
There is increasing evidence that mantle plumes should be thermally [20] and chemically [19] heterogeneous structures, and such heterogeneity should produce obvious variations in geophysical observations. Therefore, the most convincing evidence for the mantle plume hypothesis comes from geophysical methods. The thermal structure of a mantle plume is marked by a decrease in the seismic-wave velocity [19], [22], especially when a broad head forms beneath either the LAB or the 660-km discontinuity; depression or uplift perturbations in the phase transition interfaces of mantle minerals (410 and 660 km) are also marked [23], [24]. Thus, current geophysical investigations of mantle plumes mainly focus on low-velocity seismic anomalies, the delay times of teleseismic waves, and the topography of key mantle interfaces [22], [25], [26], [27], [28], [29], [30], [31]. Meanwhile, more properties are required for a further study of mantle plumes.
The electrical conductivity of mantle minerals is also sensitive to thermal variations caused by mantle plumes [32]. Therefore, constraints on the conductivity structure of the mantle via EM methods can be instrumental in mantle plume detection. The MT method has been applied to locate mantle plume heads and pathways [33], [34]. For example, MT imaging has detected the obviously enhanced conductivity of the mantle plume in the upper mantle beneath Yellowstone [35], [36]. However, petrological studies and seismic imaging have indicated that many mantle plumes, including the Yellowstone plume, likely originate from either the lower mantle or the large slow-shear velocity provinces near the CMB [17], [19], [37], [38], [39], [40], which are much deeper than the maximum penetration depth of MT soundings. Geomagnetic depth sounding (GDS), which measures EM waves possessing multiday periods [41], has become another reliable tool for mantle plume detection through its ability to elucidate the conductivity structure of the MTZ and lower mantle. However, GDS has rarely been used to provide evidence of mantle plumes, and the ability to reconstruct the electrical structure of mantle plumes has not been fully recognized to date. In this article, we conduct a series of GDS inversion tests based on mantle plume models, whereby we investigate the configuration of different station coverage areas, station spacings, noise levels, and response period ranges to discuss the potential of employing GDS to accurately detect and reconstruct mantle plumes.

A. Electrical Structure of Mantle Plumes
Mantle plumes can be regarded as adiabatic in an ambient mantle although a change in the temperature gradient should be considered within a given plume. We use the five typical mantle plume models (R1a, R1b, R1c, R2, and R3) presented by Maguire et al. [42] who established pure thermal plume models with various shapes that were produced by geodynamic simulations under different conditions. The diameter of the thin tail in the lower mantle is ∼400 km. The plume diameter increases after coming into contact with the MTZ, forming a broad head with a diameter of up to 1000 km. The head diameter just beneath the LAB is slightly larger, at ∼1200 km. These models simulate mantle plumes that are consistent with previous geophysical observations (mainly seismic-wave velocities), thereby highlighting the reliability and significance of the presented research in terms of imaging realistic mantle plume structures.
Under the condition that the electrical conductivity of mantle minerals is only sensitive to temperature, the conductivity is expressed as follows: where σ is the mineral conductivity, σ 0 is the preexponential factor, H is the activation energy, k is Boltzmann's constant, and T is the mineral temperature (unit: K). Based on (1), the corresponding conductivity of a mantle plume, σ P , with an elevated temperature, T + T , can be calculated as follows: where σ 1−D is the 1-D global average conductivity [43], [44] used as the background model (see Fig. 1).
The H value varies with the mineral assemblage and the properties of the different Earth layers. High-temperature and high-pressure conductivity experiments of the main mantle minerals have yielded H = 0.70 eV in the lower mantle [45], [46]; 1.49 and 1.93 eV in the upper and lower MTZ, respectively [47], [48]; and 1.96 eV in the upper mantle [48]. Therefore, the electrical conductivity of plumes can be calculated from (2). The conductivity is slightly enhanced (up to ∼2 times) in the lower mantle, considerably enhanced in the MTZ (up to ∼5-8 times), and heavily enhanced (up to >10 times) in the upper mantle based on the relationship between the mineral electrical conductivities and temperature changes in Earth's spherical layers (see Figs. 1 and 2).
Mantle plume models R1b and R3 demonstrate obvious interactions between the plumes and MTZ, with the R1b results providing the best fit to our current understanding of mantle plumes (see Fig. 2). For example, a mantle plume originates in the deep lower mantle with a narrow tail and develops a broad head beneath the 660-km interface and LAB. R1b is used as a typical plume model in the subsequent analysis.

B. GDS Inversion Theory
The induction source in GDS is magnetospheric equatorial ring currents [41]. These currents are concentric with the magnetic equator of Earth; the numerical simulations are, therefore, developed based on the geomagnetic spherical coordinate  [42] and shown on (left), and the increase in electrical conductivity (σ ) relative to the background conductivity (σ 0 ) is shown on (right). Plume models R1a, R1b, and R1c are snapshots that represent different stages under the same dynamic simulation, whereas plume models R1, R2, and R3 are run under changing simulation conditions. system. If we assume that the induction source is described by a single spherical harmonic function P 0 1 [41], [49], [50], then the C-response is usually estimated as follows: where a 0 is the average radius of the Earth, ω is the angular frequency, θ is the geomagnetic colatitude, and H r and H θ are the vertical (pointing downward to the center of the Earth) and colatitudinal (pointing to magnetic north) components of the magnetic field (H) on the surface, respectively. The induced geomagnetic signals that are recorded at the surface have a wide period range that spans from several days to >100 days. Equation (3) shows that the C-responses should be calculated from H. If we assume a positive time harmonic dependence of the form e iωt , then H obeys the following: where ρ is the reciprocal of the electrical conductivity σ , µ 0 is the vacuum magnetic permeability, and i is the imaginary unit. Equation (4) is solved using the staggered-grid finite difference method in a spherical coordinate system [51]. The model that is parameterized for the calculation includes the resistive air and the conductive Earth, with the outer boundary of the air placed 2a 0 from the surface and the resistivity of air set to 10 10 ·m. The inner boundary is the CMB because of the superconductive core [41]. The tangential components of H at the boundaries are specified such that (4) is valid throughout the entire model domain and the resultant numerical system remains acceptably well-conditioned. The P 0 1 source is located at a radial distance of 10a 0 from the Earth's surface to ensure that the secondary magnetic field induced by the conductive Earth could be considered negligible. A variant of the biconjugate gradient and an iteration method is used to obtain the discretized solution of (4). The divergence correction [44] is also applied to ensure the conservation of H during each iteration.
The GDS inversion can be generally expressed as an optimization problem where the objective function, (m, λ), is defined as d (m) and m (m) are the data misfit and model roughness, respectively; λ is a regularization parameter that represents the tradeoff between d (m) and m (m); and m is the electrical conductivity vector, which represents the conductivity in each prism in the case of a 3-D inversion [52]. The discretized conductivity in each grid can characterize a 3-D model with arbitrary conductivity anywhere.
Using the notation of the L p -norm measurement of the objective function, (6) is expressed as where d is the data vector, m 0 is the prior model, ψ is the forward mapping operator for calculating the model (m) responses, C −1 d is a diagonal matrix with data covariance as its diagonal elements, and C −1 m is a coefficient matrix that relates the conductivity of each grid to that of the adjacent grids in the X -, Y -, and Z -directions. We directly adopted a parameterization and smoothing method based on the model grid because it can show local and abrupt conductivity changes, which can be used to describe the mantle plume better.
The L 2 -norm inversion ( p = 2) is traditionally used in EM induction since it is generally effective in obtaining a model with responses that fit the observed data well. However, if an outlier is present in the data, then the L 2 -norm measurement of the data misfit term considerably enhances the contribution of this outlier to the objective function, thereby biasing the estimate [53]. The inversion will subsequently try to fit any outliers, leading to a deviation of the inverse model from the true model. The L 1 -norm inversion ( p = 1), in which the data misfit term is measured by L 1 -norm measurement, has a conspicuous advantage over the L 2 -norm in down weighting the impact of the outliers on the inverse model [54], thereby avoiding bias in the inverse model due to overfitting the outliers. Therefore, for synthetic data generated by adding moderate Gaussian noises, the data misfit term is measured by L 2 -norm measurement. Considering that many of the estimated responses at geomagnetic observatories are with outliers [55], synthetic data with outliers are inverted by L 1 -norm inversion to highlight the capability of the L 1 -norm inversion to deal with data with outliers.
The optimization of (5) in a GDS inversion is a nonlinear process. Here, we select the limited-memory quasi-Newton method (L-BFGS) to seek the objective function solution, which has been widely used in EM induction exploration [56]. L-BFGS is a modified form of the quasi-Newton method, which has the following basic iteration formula: where j is the number of iterations, α j is the searching step, and p j is the searching direction. p j is defined as follows: where B j is the Hessian matrix approximation and where the superscript T denotes the transpose operation. The Hessian matrix approximation avoids the direct calculation of the Hessian matrix, thereby tremendously reducing the computational requirements of the inversion [57]. The computation of the objective function gradient in (10), with respect to the model parameters at the jth iteration, includes two parts, ∇ d and ∇ m . ∇ m has an analytical expression that can readily be calculated as follows: The gradient of the data prediction error ∇ d can be computed via the chain rule as R in (11) and (12) is obtained as follows: The gradient computation in (12) requires the Jacobian matrix and forward responses of m j . The computation of the latter has been described previously. Nominally, the Jacobian matrix can be directly computed; however, a more feasible approach, such as the adjoint forward technique, can be considered [44], [58]. The adjoint forward technique allows us to compute the product of the matrix and data vector as a series of adjacent forward operations to greatly reduce the computational requirements.
In this study, the model domain is parameterized to a 1 • × 1 • grid to ensure numerical modeling accuracy. A heterogeneous grid is densified with a fineness of 2 • × 2 • horizontally to discretize the model near the target, and the size of the cell increases gradually outside the area during inversions. The configuration ensures the accuracy and efficiency of iterative inversion. All calculations are run in a computer with 16-GB RAM and eight-core Intel 1 Core 2 i7-10700 CPU in the Windows 10 platform, and the following inversion tests are all finished in 2 h.

III. INVERSION TESTS
The EM wave used in GDS analysis is only sensitive to the electrical structure in the area directly beneath and adjacent to the GDS stations. Therefore, the GDS stations should be arranged across the area where the research targets or mantle plumes are expected to be present. We will perform a series of inversion tests using different station configurations (including station coverage and spacing) based on plume R1b. We note that both the quality of the measured EM field and period ranges of the C-response have significant impacts on the detectability of conductivity anomalies, so their influence on plume detection will also be discussed. The synthetic data for the inversion tests are obtained by adding Gaussian random noise to the numerical modeling responses. The ocean effect on C-responses caused by the distribution of sea and land [59], [60], [61], [62], [63] is not considered in our theoretical inversion tests because the effect can be corrected by a surface layer comprised of sea and land in the inversion of field data [43], [55]. The ocean and land conductance can be obtained from a global surface conductance model [64].

A. Mantle Plume Responses
We calculate the C-response variations caused by plume R1b and investigate their effect on the plume responses. (180 • E, 40 • N) is chosen to be the center of mantle plumes in our tests under the consideration of convenience for modeling and calculation. The calculation costs approximately 2 min, and the variations are taken as the difference between the R1b and background model responses at a given period, as shown in Fig. 3. The maximum response variations occur at the plume edges in the latitudinal direction, thereby providing constraints on the plume boundary. The geometric center of the plume appears to be defined by the approximate location where the latitudinal response variation is zero. The response variations decrease with distance from the plume boundary, with a very weak difference observed at the largest period (113 days) and a distance of 10 • (170 • E and 190 • E in longitude, 30 • N and 50 • N in latitude). This result is significantly smaller than the EM response error (∼5%) that is commonly used in GDS inversions. Therefore, a 20 • × 20 • area that is centered on the plume axis is taken as the maximum area of station coverage in subsequent inversion tests.

B. Test A: Station Coverage
We can choose the appropriate GDS station coverage for detecting a mantle plume based on the petrological and geophysical results since few EM studies have detected and confirmed the existence of mantle plumes. The numerical modeling results for plume R1b indicate that the GDS station  coverage should be smaller than 20 • × 20 • . We employ station coverage areas of 20 • × 20 • , 16 • × 16 • , 10 • × 10 • , and 6 • × 6 • to investigate the optimal station coverage for plume detection. We set the station spacing to 2 • for each of the abovementioned areas and add 5% Gaussian random noise to the forward responses to obtain more realistic synthetic data for the inversion.
The root mean square (rms) errors of data misfit are shown in Fig. 4, meaning that the inversions have worked well, which can also be obtained by the following inversion tests. Inversion results for each test over the 250-900-km depth range are displayed in Fig. 5. The modeled extension in the longitudinal direction may be caused by the colatitude term (θ ) in the C-response calculation, which agrees with the C-response variation modeling results (see Fig. 3). The plume characteristics are well reproduced for each of the tested station coverage areas. The modeled electrical structures are closely related to the plume size in the 410-900-km depth range, which corresponds to the most sensitive GDS depth range. It is worth noting that the recovered conductivities are lower than that of the plume R1b. This is likely due to the limitation of EM imaging, whereby the EM induction associated with the extremely small-scale plume structure in the 410-520-km depth range is more sensitive to anomaly conductance rather than conductivity [65]. The inversion tests that best recover the plume are those with station coverage areas of 20 • × 20 • , 16 • × 16 • , and 10 • × 10 • , with a centrosymmetric shape obtained in cases. However, the large number of GDS stations in the first two tests compared with those for the smaller coverage areas indicates that these configurations are not our best choice in a practical (real-world) sense. Although the 6 • × 6 • station coverage area yields a result that is similar to that with the larger coverage areas, the larger, broadened anomalous zone highlights that this smaller station coverage area is unable to effectively image the plume outside of the station coverage area, with the surrounding conductive anomalies having a greater impact on the model results. Therefore, we select the 10 • × 10 • station coverage area as the best choice since it provides a good balance between using a more practical GDS station distribution and capturing the plume structure from the lower mantle to the upper mantle.
The obtained conductivities in GDS inversions are slightly lower than the real conductivity of the plume for a given layer. This conductivity leakage is largely related to the small plume size, especially in the upper and lower mantles, which makes distinguishing this small-scale anomaly from the largewavelength EM waves difficult. This reduction is essentially an inherent feature of EM induction, thereby resulting in a relatively weak vertical resolution at greater depths. The small size and electrical conductivity variations of the plume tail, in combination with its position in the lower mantle, which exceeds the GDS resolution depth, therefore, lead to poor recovery of the plume tail.

C. Test B: Station Spacing
Reducing the station spacing implies that more stations will be available to acquire data containing information on the structures of the Earth's interior and subsequently be included in the inversion, thereby improving the structural resolution [44] and more accurately revealing the mantle structure. The station spacing, therefore, plays an important role in achieving the inversion results.
Synthetic data are generated in a similar way to that outlined in Test A. Here, we employ the same GDS station coverage area (10 • × 10 • ) according to the conclusion derived from Test A and perform inversions (see Fig. 6) with station spacings of 1 • , 2 • , 5 • , and 10 • . In order to demonstrate the suitable configuration of station spacing for reconstructing the mantle plume, tests on station spacings with a 20 • × 20 • coverage (see Fig. 7) are also conducted.
The plume can be recognized by the electrical structure near the MTZ with a slightly irregular and deformed shape at 5 • station spacing with a 10 • ×10 • coverage area. The shape of the reconstructed plume has more obvious deformation, and the amplitude of conductivity is much smaller for the 10 • station spacing, especially for the 20 • × 20 • coverage area as the distribution of stations is unreasonable compared to that of 10 • × 10 • coverage area. Therefore, reducing the total number of stations by increasing the station spacing leads to less information on the plume being supplied to the inversion and weakened constraints on the inverted plume model, resulting in poor plume recovery. The accuracy of the inversion results is significantly improved when the station spacing is decreased. The inversion results for the 1 • , 2 • , and 5 • station spacings capture obvious plume characteristics and yield similar results. However, we do not recommend employing a 1 • station spacing due to the requirement of large amounts of additional stations compared with the 2 • and 5 • station spacings. The deformation in the shape of the plume [see Fig. 6(c)] presented in a 10 • × 10 • area with  Fig. 7(c)]. Therefore, either the station spacing reducing as much as possible over a smaller station coverage area (such as 2 • station spacing with 10 • × 10 • coverage area) or a larger station spacing over a larger coverage area (such as 5 • station spacing with 20 • × 20 • coverage area) should be employed during GDS data acquisition to improve mantle plume recovery. However, the GDS method is limited to the mid-latitude zone because the geomagnetic field is strongly affected by the auroral currents in high latitudes and equatorial electrojets in low latitudes [66]. Therefore, smaller station coverage is more suitable for plume detection. Fig. 3 shows that the maximum C-response variations for plume R1b can approach >10% at smaller periods (∼10 days).

D. Test C: Noise Level
These variations rapidly decrease away from the edges of the mantle plume, thereby indicating that a high signal-to-noise ratio is required for the GDS inversion. However, the observed GDS data are often contaminated with complex instrumental and environmental noise sources, such as magnetic storms, tides, human activities, heterogeneities in the lithosphere and shell, and the local ocean and continent distribution [59]. It is, therefore, difficult to ensure good data quality, and a relatively low signal-to-noise ratio of the responses is commonly present in the obtained data [55]. Here, we produce synthetic data with different noise levels to investigate the influence of noise on the inversion results.
The same station coverage area (10 • × 10 • ) and station spacing (2 • ) conditions are employed for each test, with synthetic data obtained by adding 1%, 3%, 5%, and 8% Gaussian random noise to the modeling responses of plume R1b to represent excellent-, good-, general-, and poor-quality data, respectively (see Fig. 8). The inversion results in Fig. 9 Fig. 8. Variations in the theoretical random noise data that were generated with different noise levels. The Gaussian noise level of (a)-(d) is 1%, 3%, 5%, and 8%, respectively.
indicate that the synthetic data with 1%, 3%, and 5% noise levels can all recover the mantle plume well. The inversion results for the synthetic data with 8% noise are distorted by obvious deformations and false anomalies, thereby making it difficult to restore the true shape and thermal state of the mantle plume. The results indicate that, when the noise levels either approach or exceed the conductivity variation generated by the plume itself, the signal containing the plume information becomes overprinted by these large noise levels to the point that the inversion cannot extract useful information from the noise. However, the plume shape is largely restored, and the false anomalies contained in the noise are suppressed when the L 1 -norm inversion, in which data misfit is measured by L 1 -norm instead of L 2 -norm, is adopted for these noisier data [see Fig. 9(e)]. Therefore, the quality of the GDS data should be a top priority during the data acquisition process to maximize the potential for obtaining reliable plume properties. The L 1 -norm inversion should be considered when poorquality GDS data are acquired, as this inversion approach greatly reduces the impact of large noise levels on the inversion results.

E. Test D: Period Range
The EM detection depth depends on the period range of the C-responses. Although a longer observation time will improve both the stability and maximum period of the C-responses, it comes at the expense of additional observation times and acquisition costs to obtain longer periods of the responses. Constraints on the period range that is suitable for plume detection are, therefore, strongly desired.
The synthetic data are obtained by adding 5% noise to the modeling responses. The inversion results indicate that the model from the C-responses with a period up to ∼20 days [see Fig. 10(c)] almost coincide with those from a period up to ∼100 days [see Fig. 10(a)] above the lower mantle. In the lower mantle, the reconstruction of the plume tail is limited to the topmost lower mantle, and the recovering capacity stays even in increasing periods from ∼40 to ∼100 days. This is mainly caused by the narrow diameter and weak conductance of the plume tail, and the depth of the tail in the deep lower mantle has exceeded the investigation depth of GDS. The test indicates that a period of ∼20 days is sufficient to detect the structure of the plume in the MTZ, and only the portion of the plume in the topmost lower mantle can be resolved additionally even for larger-period responses. Therefore, the effective GDS detection depth for a mantle plume is near the MTZ (∼250-660 km depth). This allows the electrical structure of the mantle plume to be effectively detected using a relatively small period range, possibly no more than up to 20 days.
The C-responses with a maximum period of ten days are sufficient to capture the basic shape of the plume when the data quality is good (e.g., the noise level is 3% or better) [see Fig. 11(d)]. However, the small maximum period means that the electrical structures in the lower mantle are poorly imaged. The overall structure is well obtained by the responses with larger periods, whereby additional information, some of which is comparable to the noise level, is included. The inversion results for a maximum period of ten days are deformed due to the small maximum period and higher noise levels, thereby highlighting that insufficient data are available to constrain the mantle plume. Furthermore, the inversion results are limited to the upper mantle and MTZ. These make it difficult to achieve an accurate judgment on the location of the root of mantle plumes. The inversion results for the responses with a period of up to ∼20 days are basically the same as those with larger periods, which should better reflect the characteristics of the plume in Earth's layers. We, therefore, suggest that the response period of the GDS observations should be >20 days taking the inherent difficulties in controlling the noise levels and GDS detection depth into consideration.

IV. DISCUSSION
The inversion tests have shown that synthetic response data with 5% Gaussian random noise from a 10 • × 10 • station coverage area with a 2 • station spacing can effectively reconstruct plume R1b. We will run additional inversion tests for the five plume types and different mantle plume sources to determine the feasibility of employing GDS to accurately characterize these plumes.

A. Indistinguishable Plume Shapes
Plume generation zone heterogeneity, the LAB, the 410-km and 660-km phase transformation discontinuities, and the plume development stage will lead to different plume states and shapes. Here, we investigate the effectiveness of GDS in detecting different plume types based on the five models shown in Fig. 2.
The inversion results in Fig. 12 reveal that current GDS technology cannot detect the narrow tail in the lower mantle and the broad head just beneath the LAB. Plume R1a, which is located in the deep lower mantle, is below the GDS resolution and is not detected in the inversion results. Obvious plume features are imaged in the mantle for the other models, but the real conductivity anomalies are weakened and spread (both laterally and vertically) in the inverted models. Fig. 12 also presents that the true plume anomalies in the 100-250-km depth range are easily transferred to the 250-410-km depth range in the inversion results. This observation is due to both the response periods, which makes the GDS observations insensitive to the depth range, and the leakage of the high-conductivity to the surrounding areas. Therefore, attention should be paid to the overall conductivity structures when interpreting the GDS inversion results, as the detailed structures of the anomalies in the shallower and deeper mantle may not reflect the true structure of the plume.
The recovered areas with increased conductivity for varied mantle plumes are mainly concentrated near the MTZ, with an overall similar structure. The range and amplitude of the high-conductivity anomalies are strongly related to the plume diameter. The plumes recovered well in the lower MTZ (520-660 km) and the topmost lower mantle (660-900 km), with a larger scale and higher conductivity revealed in the inversion results. We, therefore, suggest that the GDS inversion results can be used to assess the basic structure of the mantle plume, and the lower MTZ and the topmost lower mantle should be the best layers for discussing the property of plumes.
The conductivity values obtained by the GDS inversion are weaker than those of the real plume. The recovery rates in Fig. 11 demonstrate that the inversion results are reliable, with the recovery rate of the average conductivity in the MTZ generally exceeding 50% and that in the topmost lower mantle exceeding 80%. The maximum conductivity is slightly lower than that of the plume. Therefore, the inverted average conductivity is more suitable for constraining the range and nature of the true anomalies, while the maximum conductivity can be used to define the center and upper limit of the electrical structure, as EM induction is more sensitive to the conductance (conductivity × volume).

B. Irresolvable Plume Root
Mantle plumes generally originate from the CMB, but upwellings from the MTZ and upper mantle are also identified [19]. Therefore, we ran inversion tests for mantle plume sources at the 660-and 410-km discontinuities to investigate the feasibility of using GDS data to distinguish different mantle plume sources. We simply removed the conductivity anomalies below 660 km and 410 km in plume R1b to obtain new upwelling models, which are denoted P660 and P410, respectively. We subsequently added 5% Gaussian random noise to the forward responses and ran the inversions, with the results presented in Fig. 13.
The inversion results possess relatively obvious differences among the three mantle plume source scenarios. The anomalies obtained from the P660 and P410 inversions [see Fig. 13(c) and (d)] are concentrated above 660 km depth, and the difference between the two results in the conductivity of the lower MTZ leads to a potential but tough identification of the plume source. Simultaneously, the reconstructed anomaly in the lower mantle can be confirmed to be caused by a plume with a lower mantle origin. In brief, only the mantle plume originating from the lower mantle can be convincingly recognized.
The inversion results in Fig. 13 indicate that GDS is insensitive to the electrical structure above 250 km depth. Therefore, MTZ and the topmost lower mantle play a key role in detecting plumes, which is associated with the sensitive depth range of GDS [44], [55]. The determination of the root of plumes is arduous work for GDS, and it should be a comprehensive consideration with other geophysical and geochemical observations.

C. Influence of Record Length on Detection Depth
The GDS detection depth is determined by the period range of the C-responses, and the estimated C-responses in the frequency domain are greatly affected by the GDS record length. We estimate the C-responses from the GDS time series using a self-reference method based on Bounded Influence Remote Reference Processing (BIRRP), as in previous studies [55], [67]. The C-responses, with 16 periods that are logarithmically distributed from 3.5 to 113.8 days, can be obtained by data processing. The Fürstenfeldbrück geomagnetic station (FUR) has a stable and persistent data record (>60 years) and is located far from the ocean such that the ocean effect on FUR is negligible [59], [60]. We analyze the FUR record to determine the period range of the C-responses that can be obtained from GDS data with different record lengths (see Fig. 13).
The oscillation and squared coherency (quality indicators of the C-response) of the responses, which are presented in Fig. 13 [59], indicate that the response stability increases with record length. A response period of up to ∼100 days can be  obtained with a stable record length of >10 years, whereas a two-year record length can be converted to stable responses with a period of >20 days, and the maximum response period is only ∼10 days for a one-year record length (Fig. 14). However, shorter record lengths possess an extremely low coherency and cannot be processed by BIRRP. This may be attributed to the P 0 1 assumption, which ignores the currents corresponding to the harmonic coefficient, with the high orders and degrees corresponding to the geomagnetic fields and responses in smaller periods (e.g., auroral currents); further details on these effects can be found in [68] and [69].

D. Plume Revealed by Field Geomagnetic Data
We collected geomagnetic data recorded in China [see Fig. 15(a)], and most of the record length lasts longer than ten years. The average station spacing is approximately 5 • with a 20 • × 20 • station coverage area. Therefore, the dataset is suitable to reveal the mantle plume as we have proved. The L 1 -norm inversion is applied to convert the estimated C-responses due to the poor quality at some stations [67], [70]. The background model shown in Fig. 1 is used as the initial model in inversion. A heterogeneous grid is densified with  To eliminate the ocean effect, the surface layer with a thickness of 12.65 km is set to describe the sea and continent [43], [55], [59], [60].
A conductivity profile of inversion results is shown in Fig. 15(b). The detected continuous conductor from the lower mantle to the upper mantle beneath South China coincides with the characteristic of a reconstructed mantle plume. The most conductive zone locates near the 660 km discontinuity, corresponding to the most sensitive depth of GDS to a mantle plume. Therefore, we speculate that the conductor should be a mantle plume originating from the lower mantle. The speculation is also proved by the seismic tomography [31], [71], [72], [73] and perturbations of 410-and 660-km discontinuities [74].

V. CONCLUSION
In this article, the optimal GDS acquisition parameters for mantle plume construction have been determined based on the inferred electrical structures of different mantle plume types, which are based on geodynamic simulations under different conditions and laboratory conductivity measurements of different minerals under the temperature-pressure conditions of the deep mantle. The inversion tests of synthetic data and field data were conducted to assess the effectiveness of GDS in reconstructing mantle plumes yielded the following key results.
1) The GDS inversion is mainly sensitive to the mantle plume structures near the MTZ, with high recovery of the mantle plume structure in this depth range and the ability to distinguish 10 • diameter plumes.
2) Increasing both the station coverage area and station density can lead to a more accurate reconstruction of mantle plumes. However, a 10 • × 10 • station coverage area with a ∼2 • station spacing is sufficient for imaging an ∼10 • diameter plume with a normal noise level in GDS data (∼5% noise).
3) The GDS inversion results suffered when large noise levels were introduced to the input data. When responses at observations contain little noise (<3% noise), a continuous one-year record of the geomagnetic field is effective to construct the mantle plume. For responses with normal noise (∼5% noise), a continuous two-year record is required. 4) Noises in larger levels (∼8% noise or more) will overprint the plume information, resulting in abnormal deformations or false anomalies in inversion results. Application of the L 1 -norm inversion, whereby the data misfit is measured by L 1 -norm, can improve the reliability of the inversion results when larger noise levels are present in the GDS observations. 5) A mantle plume originating from the lower mantle is discovered by GDS beneath South China.