A Mantle Plume Beneath South China Revealed by Electrical Conductivity Obtained from Three-Dimensional Inversion of Geomagnetic Data

A three-dimensional electrical conductivity model of the mantle beneath South China is presented using the geomagnetic depth sounding method in this paper. The data misfit term in the inversion function is measured by the L1-norm to suppress the instability caused by large noises contained in the observed data. To properly correct the ocean effect in responses at coastal observatories, a high-resolution (1° × 1°) heterogeneous and fixed shell is included in inversion. The most striking feature of the obtained model is a continuous high-conductivity anomaly that is centered on ~(112° E, 27° N) in the mantle. The average conductivity of the anomaly appears to be two to four times higher than that of the global average models at the most sensitive depths (410–900 km) of geomagnetic depth sounding. Further analysis combining laboratory-measured conductivity models with the observed conductivity model shows that the anomaly implies excess temperature in the mantle. This suggests the existence of a mantle plume, corresponding to the Hainan plume, that originates in the lower mantle, passes through the mantle transition zone, and enters the upper mantle. Our electrical conductivity model provides convincing evidence for the mantle plume beneath South China.


Introduction
Mantle plumes can transport material and energy to the surface. They may be closely related to the breakup of the paleo-continents, the formation of large igneous provinces, the eruption of a series of intraplate volcanoes, and even mass extinctions [1,2]. Therefore, the reconstruction of mantle plumes plays a key role in understanding mantle convection and Earth's evolution. The temperature and mineral compositions of a plume are obviously different from those of the surrounding mantle, and a mantle plume interacts with the mantle transition zone (MTZ) and other interfaces in Earth during its rise, causing significant changes in the structure of the mantle [3]. Therefore, the most convincing evidence for the existence of a mantle plume comes from geophysical observations of the mantle.
Geophysical methods, particularly seismology and seismic tomography, are used to detect plumes. Many achievements have been made in plume research through seismic imaging [4][5][6][7][8]. Recently, a low-velocity structure was discovered in the mantle below Hainan Island [9][10][11][12][13][14]; this structure originates in the lower mantle, crosses the MTZ, and reaches the lithosphere. This characteristic suggests that the low-velocity structure is a mantle plume. Therefore, the Hainan plume has been put forward to explain the formation and evolution of the South China Sea, the mechanism of Hainan's volcanoes, and other scientific issues. However, there are still large disputes among different research groups about the location, shape, temperature, and other properties of this mantle plume. For

GDS Theory
The inducing source for GDS is the slowly changing ring currents in the magnetosphere [17]. These currents are concentric with the magnetic equator of Earth, and so the numerical simulation is developed based on the geomagnetic spherical coordinate system. The widely used C-response of GDS is estimated from the H r (the vertical component pointing downwards towards the center of Earth) and H θ components (the colatitudinal component pointing towards magnetic north) of the magnetic field (H) at the surface. With the assumption that the ring currents can be described by the spherical harmonic function P 0 1 [17][18][19], the C-response can be calculated by where a 0 is the average radius of Earth, ω is the angular frequency, and θ represents the geomagnetic colatitude (0-180 • ). The induced geomagnetic signal collected at the surface for GDS has a period of several days to more than 100 days. Equation (1) shows that C-responses should be calculated from the magnetic field H. Under the assumption of positive time harmonic dependence of the form e iωt , H obeys where ρ is the reciprocal of electrical conductivity σ, µ 0 is the vacuum magnetic permeability, and i is the imaginary unit. Equation (2) can be solved by means of the staggered-grid finite difference method in a spherical coordinate system [20]. The model parameterized for calculation includes resistive air and conductive Earth. The outer boundary of air is 2a 0 from the surface, and its resistivity is set to a moderately large finite value of 10 10 Ω·m. The inner boundary of Earth is the core-mantle boundary (CMB) due to the superconductive core [17]. The tangential components of H at the boundaries are specified so that Equation (2) holds throughout the space domain, while the resultant numerical system remains acceptably well-conditioned. The location of the P 0 1 source is placed at a radial distance from Earth's surface of 10a 0 to ensure that the secondary magnetic field induced by the conductive Earth can be considered negligible. A variant of the biconjugate gradient and an iteration method are used to obtain the solution of discretized Equation (2) [21]. To ensure that H is conservative during iteration, the divergence correction [22] is also applied.

L-BFGS Inversion
GDS inversion can generally be expressed as an optimization problem where the penalty function Φ(m, λ) is defined by where Φ d (m) and Φ m (m) are data misfit and model roughness, respectively, and λ is the regularization parameter, which is used to trade off Φ d (m) and Φ m (m). m is the electrical conductivity vector, and in the case of three-dimensional (3D) inversion, it represents the conductivity in each prism [5].
Using the notation of L p -norm measurement of the objective function, Equation (4) is where d is the data vector, m 0 is the prior model, ψ represents the forward mapping operator used to calculate the responses of model m, W d is a diagonal matrix with data covariance as diagonal elements, and W m is a smoothing matrix used to relate the conductivity of each grid to that of the adjacent grids in three directions (X, Y, and Z directions). The correlation among adjacent cells increases as the value of the smoothing coefficients increases from 0 to 1. L p -norm inversion can be realized by assigning different values of p.
In the traditional inversion approach, differentiating both sides of Equation (5) with respect to the space domain model parameters and neglecting the higher order terms of the Taylor-series expansion allow the linear system of equations to be solved at each iteration, as follows: where where ε is a small number to ensure a solution when x = 0 and p corresponds to the L p -norm inversion. The optimization of Equation (3) in the case of GDS inversion is nonlinear, and we select the limited-memory quasi-Newton method (L-BFGS), which has been widely used in electromagnetic induction exploration [23] to seek the solution of the penalty function. L-BFGS is a modified form of the quasi-Newton method. The basic iteration formula of L-BFGS is where and where k is the number of iterations, α k is the searching step, p k is the searching direction, and B k is the approximation of the Hessian matrix [24]. The approximation of the Hessian matrix avoids calculating the Hessian matrix directly, thus tremendously reducing the requirements for computer storage and computation time.
The computation of Equation (6) requires the calculation of the Jacobian matrix and forward responses. The latter can be easily accomplished. Nominally, the Jacobian matrix can be directly computed, but a more feasible method, the adjoint forward technique, can be considered [21,22,25]. Using this method, we compute the product of the matrix and the data vector, which is split into a few adjacent forward operations, thus greatly reducing computational requirements.

Data Processing
C-responses can be estimated using the BIRRP software package [26], which is a combination of the standard M estimation method and hat matrix diagonal statistical analysis; thus, it can eliminate the interference in H r and H θ and correlated noise in both. The application of the remote reference method is suggested in BIRRP, but the self-reference method has been proven to show only a negligible difference from the remote reference method [27] in the long periods used for the present work; therefore, in this paper, we estimate C-responses by means of the self-reference method based on BIRRP.
Hourly mean value time series of the three components of the geomagnetic field can be obtained from the World Data Center. In this paper, the selected time series are from 6 years up to approximately 60 years. In addition, the Geomagnetic Network of China provides data from deployed geomagnetic stations in China, which significantly enhances the number of stations that can be used in GDS inversion. After careful selection according to the duration and noise that we presented previously [21], we obtained 15 geomagnetic observatories in South China for further consideration ( Figure 1). The details of the observatories are listed in Table 1. The estimated C-responses and their errors within 16 periods from 3.5 to 113 days at the selected stations are presented in Figure 2. The estimation of C-responses in the frequency domain strongly depends on the data length in the time domain. However, the data length varies at different stations, besides the observed data being strongly affected by environmental noise and the ocean effect, leading to a huge variability both in the real and imaginary parts of C-responses, especially for responses in a long period. The variation in responses at different stations indicates that the electrical structure beneath South China is heterogeneous. The squared coherency of C-responses, which is commonly treated as a quality indicator of C-responses, is depicted in Figure 2c, showing that most of the C-responses are of good quality, while some stations produce C-responses of poor quality. Considering the strong fluctuation for some stations, the L 1 -norm measurement is used to measure the data misfit to suppress the influence of data noise on the inversion results [21].

Influence of the Ocean
The C-responses at coastal observatories are significantly influenced by ocean induction effects (OIEs) due to the large contrast in conductivity between oceans and continents, especially for responses over short periods [27][28][29][30]. Some of our selected observatories are located near the coastline, so we examine the OIEs on C-responses. Based on the global average one-dimensional (1D) model ( Figure 3) derived from vector magnetic field data observed by satellite [31], we calculated the C-response at the Guangzhou (GZH) and Sheshan (SSH) stations near the coastline. The C-responses of the 1D model covered by a shell comprising oceans and continents at the two stations were also calculated. As shown in Figure 4, both the real and imaginary components of the C-response are significantly affected in the examined periods. The influence on the imaginary components is more obvious, and the OIEs are stronger in GZH than in SSH. The influence of the OIEs at the two stations can reach a period of 30 days if 5% variations are taken as a reference, which are commonly used as an error floor in electromagnetic induction inversion [32][33][34]. According to the relationship between electrical conductivity and the C-response, the conductivity is proportional to the imaginary component of the C-response, and the penetration depth is directly proportional to the real component. The strong variations in C-responses correspond to the extremely enhanced conductivity in the shallow mantle, which do not coincide with the actual situation. Therefore, the response of stations affected by OIEs must be corrected. mula is as follows: where is the corrected response, is the observed response at the station, and is the correction coefficient, which can be calculated from the response of the 1D model and the response of the 1D model with a covering shell comprising oceans and continents. Obviously, the correction coefficient changes with the 1D model. This makes it difficult to identify a 1D model that can be suitable for all stations because the electrical structure of Earth is unknown and three-dimensional. Therefore, in our 3D inversion of GDS, the conductivity of the shell with ocean and land was treated as the surface layer of Earth. The shell was considered in the forward numerical modeling but was fixed throughout the inversion [21,36].

Electrical Conductivity Model
The C-responses ( Figure 2) at the 15 stations in South China were inverted via 3D GDS inversion, in which the data misfit and model roughness were measured using L1and L2-norm measurements, respectively. To reduce the influence of the electrical structures of the surrounding mantle, stations near the research area were also considered in our inversion. The data error used to normalize the data misfit in inversion was obtained from the estimated responses. Data with large uncertainties should be excluded from the traditional L2-norm inversion, but in L1-norm inversion, their contributions could be suppressed by data error normalization and L1-norm measurement. Therefore, it is expected that we could obtain a reliable 3D electrical structure of the mantle beneath South China.
The initial model of our 3D L1-norm inversion adopted the 1D model in Figure 3. To match the major mineral phase transitions in the mantle, the conductivity was allowed to jump at 410 km, 520 km, and 670 km. A heterogeneous grid was densified with a fineness of 3° × 3° horizontally to discretize the model in the research area, and the size of the cell increased gradually outside the area (as shown in Figure 5). To eliminate the OIEs, a lat- The ratio method can be used for OIE correction [35], and the specific correction formula is as follows: C corr = C obs ·k, and k = C 1D /C 1D+shell (11) where C corr is the corrected response, C obs is the observed response at the station, and k is the correction coefficient, which can be calculated from the response of the 1D model and the response of the 1D model with a covering shell comprising oceans and continents. Obviously, the correction coefficient changes with the 1D model. This makes it difficult to identify a 1D model that can be suitable for all stations because the electrical structure of Earth is unknown and three-dimensional. Therefore, in our 3D inversion of GDS, the conductivity of the shell with ocean and land was treated as the surface layer of Earth. The shell was considered in the forward numerical modeling but was fixed throughout the inversion [21,36].

Electrical Conductivity Model
The C-responses ( Figure 2) at the 15 stations in South China were inverted via 3D GDS inversion, in which the data misfit and model roughness were measured using L 1 -and L 2 -norm measurements, respectively. To reduce the influence of the electrical structures of the surrounding mantle, stations near the research area were also considered in our inversion. The data error used to normalize the data misfit in inversion was obtained from the estimated responses. Data with large uncertainties should be excluded from the traditional L2-norm inversion, but in L 1 -norm inversion, their contributions could be suppressed by data error normalization and L 1 -norm measurement. Therefore, it is expected that we could obtain a reliable 3D electrical structure of the mantle beneath South China.
The initial model of our 3D L 1 -norm inversion adopted the 1D model in Figure 3. To match the major mineral phase transitions in the mantle, the conductivity was allowed to jump at 410 km, 520 km, and 670 km. A heterogeneous grid was densified with a fineness of 3 • × 3 • horizontally to discretize the model in the research area, and the size of the cell increased gradually outside the area (as shown in Figure 5). To eliminate the OIEs, a lateral grid of 1 • × 1 • of the surface layer with a thickness of 12.65 km was considered to more precisely describe the distribution of the sea and continent to generate the shell in 3 • × 3 • for calculation and ensure OIEs were accounted for with sufficient accuracy [28].

Electrical Conductivity Model
The C-responses (Figure 2) at the 15 stations in South China were inverted via 3D GDS inversion, in which the data misfit and model roughness were measured using L1and L2-norm measurements, respectively. To reduce the influence of the electrical structures of the surrounding mantle, stations near the research area were also considered in our inversion. The data error used to normalize the data misfit in inversion was obtained from the estimated responses. Data with large uncertainties should be excluded from the traditional L2-norm inversion, but in L1-norm inversion, their contributions could be suppressed by data error normalization and L1-norm measurement. Therefore, it is expected that we could obtain a reliable 3D electrical structure of the mantle beneath South China.
The initial model of our 3D L1-norm inversion adopted the 1D model in Figure 3. To match the major mineral phase transitions in the mantle, the conductivity was allowed to jump at 410 km, 520 km, and 670 km. A heterogeneous grid was densified with a fineness of 3° × 3° horizontally to discretize the model in the research area, and the size of the cell increased gradually outside the area (as shown in Figure 5). To eliminate the OIEs, a lateral grid of 1° × 1° of the surface layer with a thickness of 12.65 km was considered to more precisely describe the distribution of the sea and continent to generate the shell in 3° × 3° for calculation and ensure OIEs were accounted for with sufficient accuracy [28].  The 3D inversion started with an initial regularization parameter of 1.0. During iterations when the relative change in the data fitting error was less than the given threshold value, the regularization parameter was halved. The inversion process continued iteratively until the target data misfit was reached or the regularization parameter was rather small. The inversion was repeated for a number of different regularization parameters to better assess the stability of any detected mantle anomalies. The results of these inversions show that the distribution of anomalies is similar for different inverse models, except that the boundary sharpness is different. After 59 iterations, the inversion terminated with a regularization parameter smaller than 10 −4 . The root mean square (RMS) of data misfit is 1.81, which is larger than the expected value of 1.0. This large RMS is due to some unreliable data included in the inversion. If we exclude the responses whose RMS is larger than 5.0, the RMS decreases to 1.19. When using L 1 -norm inversion, the models from the two inversions share almost the same electrical structures. The results of inversions with an initial regularization parameter equaling 100 and stations excluding responses with a large RMS are shown in Figure 6 to strengthen the reliability of our preferred model. The data the boundary sharpness is different. After 59 iterations, the inversion terminated with a regularization parameter smaller than 10 −4 . The root mean square (RMS) of data misfit is 1.81, which is larger than the expected value of 1.0. This large RMS is due to some unreliable data included in the inversion. If we exclude the responses whose RMS is larger than 5.0, the RMS decreases to 1.19. When using L1-norm inversion, the models from the two inversions share almost the same electrical structures. The results of inversions with an initial regularization parameter equaling 100 and stations excluding responses with a large RMS are shown in Figure 6 to strengthen the reliability of our preferred model. The data fitting curves are plotted in Figure 7. The curves show a good fitting for most responses, implying a convincing result for our inversion.      (Anomaly A) is the enhanced conductivity area, which is centered on( 112 • E, 27 • N). The anomaly is nearly vertical and continuous. It can be found in the lower mantle, MTZ, and upper mantle. The strongest variation in conductivity occurs in the lower MTZ, with conductivity reaching approximately 7 S/m at the center. However, this extremely high conductivity is not suitable for further analysis, because the electromagnetic induction is primarily sensitive to the integrated conductance of a conductivity body which prefers to generate a compensating higher conductivity in the core [37]. In the topmost lower mantle (670-900 km), the anomaly is similar to that in the lower MTZ both in shape and conductivity value. The average conductivity of Anomaly A appears to be two to four times higher than that of the global average models at the most sensitive depths (410-900 km) of GDS. The variation in conductivity is much weaker in the upper MTZ (410-520 km), but enhanced conductivity is observed in the shallow upper mantle. Although Anomaly A seems to extend to the deeper lower mantle (900-1200 km), we could not confirm this extent because these depths were outside of the sensitive zone of GDS, and the anomaly may have been caused by leakage from the conductive zone at 670-900 km [22]. Therefore, the anomaly at these depths will not be discussed here.

Robustness of Anomaly A
To determine whether Anomaly A is strictly required by the data, we repeated the L1-norm inversion with the same control parameters but set a hard prior bound in the region where the cell conductivity is fixed to that of the starting model and not permitted to vary during the iterative inversion as a free parameter. The results of the re-inverted test are shown in Figures 9 and 10. Figure 9 shows that if the area of Anomaly A is restricted, the inversion produces new conductive areas around the periphery of the previously conductive zone (rectangle with red dashed line), and the new areas are the most obvious at depths of 670-900 km.

Robustness of Anomaly A
To determine whether Anomaly A is strictly required by the data, we repeated the L 1 -norm inversion with the same control parameters but set a hard prior bound in the region where the cell conductivity is fixed to that of the starting model and not permitted

Discussion
The conductivity of the mantle is affected by the mineral composition, temperature, and volatiles such as water [38,39]. The conductivity values of minerals measured under high-temperature and high-pressure conditions with varying water contents in the laboratory [40] allow us to explore the properties of Anomaly A. Heterogeneity in composition and temperature is mainly created by subducting slabs and mantle plumes, while water is carried and released by subducting slabs. The contribution of water can be excluded because there is little evidence of a slab that penetrates or stagnates in the MTZ beneath South China. This exclusion can be strengthened by the elevated 660 km discontinuity at

Discussion
The conductivity of the mantle is affected by the mineral composition, temperature, and volatiles such as water [38,39]. The conductivity values of minerals measured under high-temperature and high-pressure conditions with varying water contents in the laboratory [40] allow us to explore the properties of Anomaly A. Heterogeneity in composition and temperature is mainly created by subducting slabs and mantle plumes, while water is carried and released by subducting slabs. The contribution of water can be excluded because there is little evidence of a slab that penetrates or stagnates in the MTZ beneath South China. This exclusion can be strengthened by the elevated 660 km discontinuity at  Figure 9 shows that if the area of Anomaly A is restricted, the inversion produces new conductive areas around the periphery of the previously conductive zone (rectangle with red dashed line), and the new areas are the most obvious at depths of 670-900 km. The conductivity of the zone from 900 to 1200 km increases to a rather large value compared to that of the preferred model. Figure 10 shows the fitting curves of C-responses for three observatories (CDP, GZH, and WHN) located near the region of Anomaly A. It is evident that when suppressing the conductive feature by setting a hard prior bound (i.e., the conductivity in this zone must be equal to the background conductivity), the model response misfit deteriorates significantly at the three stations, particularly at GZH and WHN. The response shifts considerably from the observed data at nearly all periods, meaning that the conductivity at all depths is inconsistent with the real mantle situation. We therefore conclude that Anomaly A is a required feature of the data.

Discussion
The conductivity of the mantle is affected by the mineral composition, temperature, and volatiles such as water [38,39]. The conductivity values of minerals measured under high-temperature and high-pressure conditions with varying water contents in the laboratory [40] allow us to explore the properties of Anomaly A. Heterogeneity in composition and temperature is mainly created by subducting slabs and mantle plumes, while water is carried and released by subducting slabs. The contribution of water can be excluded because there is little evidence of a slab that penetrates or stagnates in the MTZ beneath South China. This exclusion can be strengthened by the elevated 660 km discontinuity at the location measured using the receiver function technique [41], which is contradictory to the depression caused by water [42]. The influence of a plume on conductivity is more dependent on the excess temperature rather than on compositional differences. Therefore, the enhanced conductivity of Anomaly A may be a result of high temperatures.
In the lower mantle, the conductivity of the lower mantle can be estimated by where σ is the electrical conductivity of the lower mantle, σ 0 is the pre-exponential factor, and k is the Boltzmann constant. P and T are the pressure and temperature, and their values at different depths in the mantle can be extracted from Xu et al. [43]. ∆E and ∆V are the activation energy and activation volume, respectively. Taken from the literature [43], σ 0 = 74 S/m, ∆E = 0.7 ev, and ∆V = −0.55 ± 0.01 cm 3 /mol. In the uppermost lower mantle, the average temperature is approximately 1900 K, which can be calculated from the global average conductivity according to Equation (12) and matches the result presented by Xu et al. [43]. To fit the average conductivity of Anomaly A, a temperature at approximately 2300 K, which is 400 K higher than the global average, is required ( Figure 11b). Additionally, the temperature at the CMB is approximately 2550 K; assuming an adiabatic upwelling mantle plume originating from the CMB, this estimate is reasonable. In the MTZ, wadsleyite and ringwoodite have a relatively high water capacity [44,45]. We can estimate the average water contents in the MTZ according to the global average 1D model and the electrical conductivity model of Yoshino et al. [38] where C w is the water content, σ 0 is the pre-exponential factor, α is a constant accounting for geometrical factors, and H is the activation enthalpy. Subscripts H and P denote small polaron and proton conduction, respectively. The values of parameters contained in Equation (13) measured by Yoshino et al. [38] are shown in Table 2. The estimated average water content in the MTZ is used to evaluate the temperature of Anomaly A. Figure 11b shows that in the lower MTZ, the conductivity of ringwoodite with a temperature approximately 300 K higher than the geotherm is close to the average conductivity of Anomaly A. In the upper MTZ, a comparable temperature approximately 200 K higher than the geotherm is obtained based on the average conductivity.
[43], = 74 S/m , Δ = 0.7 ev , and Δ = −0.55 ± 0.01 cm /mol. In the uppermost lower mantle, the average temperature is approximately 1900 K, which can be calculated from the global average conductivity according to Equation (12) and matches the result presented by Xu et al. [43]. To fit the average conductivity of Anomaly A, a temperature at approximately 2300 K, which is 400 K higher than the global average, is required ( Figure  11b). Additionally, the temperature at the CMB is approximately 2550 K; assuming an adiabatic upwelling mantle plume originating from the CMB, this estimate is reasonable. Figure 11. The conductivity cross-section along 27° N (a) and the rock physics model used to interpret the conductivity of Anomaly A (b). The conductivity of Anomaly A is delimited in the light red zone, with the red line denoting its average conductivity. The black lines with labeled temperatures represent the conductivity of minerals with corresponding temperatures.
In the MTZ, wadsleyite and ringwoodite have a relatively high water capacity [44,45]. We can estimate the average water contents in the MTZ according to the global average 1D model and the electrical conductivity model of Yoshino et al. [38] where is the water content, is the pre-exponential factor, is a constant accounting for geometrical factors, and is the activation enthalpy. Subscripts H and P denote small polaron and proton conduction, respectively. The values of parameters contained in Equation (13) measured by Yoshino et al. [38] are shown in Table 2. The estimated average water content in the MTZ is used to evaluate the temperature of Anomaly A. Figure  11b shows that in the lower MTZ, the conductivity of ringwoodite with a temperature approximately 300 K higher than the geotherm is close to the average conductivity of Anomaly A. In the upper MTZ, a comparable temperature approximately 200 K higher than the geotherm is obtained based on the average conductivity.  In the upper mantle, the conductivity is mainly determined by olivine. The conductivity of hydrous olivine has been measured [46] and can be calculated by where σ 0 and ∆H are the pre-exponential factors and activation enthalpies, respectively. The superscripts vacancy, polaron, and hydrous represent the conductive mechanisms in olivine. Values of the parameters measured by Gardés et al. [46] are listed in Table 3. Anomaly A in the deep upper mantle can be explained by a temperature enhanced by approximately 200 K in the region. Through the above analysis, Anomaly A is seen to be caused by high temperature. The temperature is highest in the lower mantle and gradually decreases in the MTZ and upper mantle. The correlation in the location of this structure in these different zones indicates that the increased temperatures have a common origin. Therefore, we speculate the existence of a mantle plume beneath South China, which is commonly defined as the Hainan plume. As shown in our electrical structure, the Hainan plume originates from the lower mantle. The low viscosity of the high-temperature mantle plume causes it to rise rapidly in the lower mantle, forming a narrow tail that is beyond the resolution of GDS. When the plume impinges on the MTZ, the exothermic reaction caused by the phase transformation at the 670 km discontinuity interface can block upwelling and result in the accumulation of a tremendous volume of hot plume materials beneath the MTZ, forming a large head [47]. Simultaneously, the bottom of the MTZ is heated by the impeded hot plume, significantly enhancing conductivity in the lower MTZ. The broad head and heated lower MTZ are revealed in GDS images, as shown in Figures 7 and 11a. In the MTZ, the viscosity differences between the plume and wadsleyite and ringwoodite allow the plume to pass through the MTZ rapidly through a narrow channel. Therefore, it is difficult to detect the plume itself if the MTZ heated by the plume is ignored. This is the reason why only slight conductivity increases are observed in the upper MTZ (Figure 11a). The corresponding relationship between location and temperature between the upper mantle and the lower mantle is a clear indicator that the Hainan plume passes through the MTZ and enters the upper mantle. However, limited by the detection depth of GDS, we are unable to trace the rising plume in the upper mantle.
Our deduction is also supported by previous seismological images. Xia et al. [10] discovered a continuous low-velocity anomaly in and around the MTZ, which extends down to the lower mantle deeper than 1100 km, as marked by the black circles in Figure 7. Considering that the electrical structure not covered by geomagnetic stations cannot be well constrained by GDS, we only show the seismic structure in the region covered by our selected stations. The low-velocity anomaly is suggested to be the tail and head of the Hainan plume, which originates in the lower mantle and feeds the Hainan hotspot. Other wave velocity structures [48][49][50] and perturbations of 410 km and 660 km discontinuity [41] also suggest the existence of the Hainan plume. The high temperature associated with the Hainan plume can concurrently decrease the seismic velocity and increase the electrical conductivity significantly [3].

Conclusions
In this paper, we have presented the first 3D electrical structure of the mantle beneath South China from geomagnetic stations in the area. The L 1 -norm inversion method was used to suppress the influence of data with considerable noise, and model roughness was measured by an L 2 -norm as in the traditional inversions to obtain a smooth model. The self-reference method based on BIRRP was used to obtain the C-response data used in our inversion in periods from 3.5 to 113 days, making GDS sensitive to conductivity at depths of 250-1600 km. The inverted model revealed the presence of a vertical and continuous conductivity anomaly in the depth range of 250-900 km centered on (112 • E, 27 • N).
The results of the laboratory measurements of the electrical conductivity of mantle minerals at high temperatures and high pressures were used to interpret the nature of Anomaly A. This analysis strongly suggests that the anomaly is caused by the high temperature of the Hainan plume. The electrical conductivity model shows that the Hainan plume originates in the lower mantle and that it passes through the MTZ rapidly and ultimately enters the upper mantle. We provide convincing evidence for the existence of the Hainan plume based on the electrical structure of the mantle converted from the geomagnetic data in South China via 3D GDS inversion.
The detection depth of GDS limits its ability to distinguish the Hainan plume in the shallow upper mantle. Therefore, appropriate methods (e.g., the long-period magnetotelluric method) can be used to assess the electrical structure in the upper mantle. Therefore, future studies with joint inversion of GDS and long-period magnetotelluric data will strongly advance the exploration of the Hainan plume.