Electrical Conductivity Evidence for the Existence of a Mantle Plume Beneath Tarim Basin

: This paper proposes using a simulated annealing (SA) calculation to perform one-dimensional inversion of Geomagnetic Depth Sounding (GDS) to obtain the conductivity information of the lower mantle beneath the Tarim area, to calculate the temperature of the lower mantle according to the relevant formula of the petrophysical experiment, and to provide evidence of the existence of the Tarim mantle plume. The data used for inversion originate from the China Geomagnetic Network Center. This article uses theoretical data to prove that the simulated annealing algorithm can invert the true conductivity model when the data do not contain noise. However, when the data contain noise, it is more accurate to use the statistical expected value of the high-quality conductivity model during the simulated annealing inversion process as the optimal conductivity model rather than the classic simulated annealing algorithm. The simulated annealing inversion results of only four stations in Tarim area show that the conductivity of the top of the lower mantle and the upper part of the mantle transition zone in Tarim area is higher than the global average, and it is speculated that the temperature is 150k–450k higher than the global average. This is important evidence for the existence of the mantle plume beneath the Tarim Basin.


Introduction
The Tarim Basin is located in the western part of China. Since 50~60 Ma, the India/Asia collision has resulted in strong uplift, shortening, and strike slip deformation around the Tarim Basin, resulting in the formation of a complex ringed mountain system [1]. The movement of the deep mantle is one of the important factors affecting the activity of the top crust. In a large number of seismic profiles, it is found that there is a large-scale transverse low velocity anomaly beneath the Tarim Basin [2,3]. Kosarev pointed out that this may indicate the existence of a mantle plume [4].
A mantle plume, first proposed by Morgan [5,6], is a columnar body produced in the thermal boundary layer of the mantle, driven by buoyancy anomalies and upwelling of deep mantle material [7,8]. Campbell believes that the mantle plume originated from the core-mantle boundary at 2900 km or the seismic discontinuity at about 660 km [9]. The diameter of the mantle plume can be up to hundreds of kilometers, and it has an abnormally high temperature. This localized high temperature anomaly can cause a dramatic increase in electrical conductivity in the transition zone and the lower mantle. Therefore, this renders us an opportunity to use the geomagnetic depth sounding (GDS) method to explore the possible existence of the mantle plume by studying the conductivity model of the mantle transition zone beneath Northwest China [10].
The electrical conductivity of the deep part of the earth, especially in the range of 410~1600 km, can be effectively reflected by the variation with periods of the C-response of GDS [11]. This parameter is estimated from the horizontal and vertical magnetic field data recorded by geomagnetic stations [12,13]. The intensity of the induced magnetic field is affected by the conductive convective fluid in the core and the mantle, crust, ionosphere, magnetosphere, and ocean current systems and changes over time [14,15]. Therefore, by investigating the changes in the intensity of the surface magnetic field, the internal conductivity of the earth at different depths can be obtained, from which the physical and chemical characteristics and the dynamic mechanisms in Earth's interior can be resolved [16,17].
Inversion is the key to convert C-responses to electrical conductivity profiles. For instance, Munch et al. used stochastic optimization and model exploration techniques to perform one-dimensional inversion of geomagnetic depth sounding, and combined with experimental data of temperature, pressure, and the water content, they obtained temperature information of the mantle transition zone in Europe and Africa [18]. Zhang et al. carried out a one-dimensional inversion of the smooth model based on the geomagnetic data obtained by the global smoothing constraint technique and obtained the water content distribution of the mantle transition zone in eastern China [19]. These one-dimensional inversion results perfectly show the conductivity characteristics of the deep part of the earth, and can provide estimations of temperature and the water content in transition zones. However, these study areas are far from Tarim area. The global GDS inversion may reflect the conductivity characteristics beneath Tarim area, for example, using the three-dimensional nonlinear conjugate gradient method, the first electrical conductivity image of mantle transition zone was obtained by Kelbert et al. [20], but they focused the high conductivity structure of the mantle beneath NE China, while Semenov et al. paid more attention to the electrical conductivity in the transition zone beneath Europe and Africa by using the limited memory quasi-Newton method [21,22]. Nonetheless, due to the sparse data set used, the resolution of the electrical structure beneath Tarim area is very low.
In this paper, we try to invert the GDS response in Tarim area by the non-heuristic algorithm of simulated annealing (SA) to provide a new reference for studying the deep electrical conductivity structure beneath the Tarim area. Similar to the genetic algorithm and artificial neural network [23][24][25], SA is a global optimization method that has been applied with good effect to geophysical inversion [26][27][28]. Traditional SA inversion accepts the model of the last iteration as the inversion result. However, this resultant model can vary significantly due to the nature of random inversion techniques, such as SA. Therefore, here we will conduct the statistical analysis of all effective random inversion models, and use the statistical model as the one-dimensional random inversion result of GDS data. To this end, the paper first briefly introduces the basic process of obtaining the C-response of geomagnetic depth sounding; then, the article addresses the stochastic inversion technique. After that, the above technique is applied to invert the GDS data from four geomagnetic stations (Wushi (WUS), Urumqi (WMQ), HongQian (HOQ), and WenQuan (WEQ)) in Northwest China to obtain the conductivities in and below transition zone. Finally, the model robustness followed by the model interpretation with abnormal temperature beneath Tarim area is discussed. Our research results can provide a very important scientific basis for the study of the evolutionary process of the orogenic belt around the Tarim Basin. From the perspective of Geophysics, we can see that there may be a high temperature mantle plume in the lower part of the Tarim Basin (this mantle plume has little to do with the Tarim large igneous province).

Data
We collected the data from the four stations in Northwest China from the China Geomagnetic Network Center ( Figure 1). The data information is listed in Table 1. These geomagnetic stations have just been established; the data have not yet been explored. The application of the new data will provide new insights into the deep structure of the earth.

Data
We collected the data from the four stations in Northwest China from the China Geomagnetic Network Center ( Figure 1). The data information is listed in Table 1. These geomagnetic stations have just been established; the data have not yet been explored. The application of the new data will provide new insights into the deep structure of the earth. In order to make better use of these data, it is necessary to use approximate methods to estimate the C-response, which can be written in the following form (1) In the formula, tanθ is regarded as the compensation term of the source space. a0 is the radius of the earth and its value is 6370 km. In the study of geomagnetic depth sounding (GDS), a spherical coordinate system whose origin is defined at the center of the earth is used. Among them, Hr points to the center of the earth, Hθ is the tangential component of the magnetic field, and θ is the geomagnetic colatitude. C-response establishes the relationship between surface observation data and the internal electrical structure of the earth [29,30].
All four geomagnetic stations record the changes of the geomagnetic component at the time. Figure 2 shows the long and continuous signals recorded by the four geomagnetic stations.  In order to make better use of these data, it is necessary to use approximate methods to estimate the C-response, which can be written in the following form In the formula, tanθ is regarded as the compensation term of the source space. a 0 is the radius of the earth and its value is 6370 km. In the study of geomagnetic depth sounding (GDS), a spherical coordinate system whose origin is defined at the center of the earth is used. Among them, H r points to the center of the earth, H θ is the tangential component of the magnetic field, and θ is the geomagnetic colatitude. C-response establishes the relationship between surface observation data and the internal electrical structure of the earth [29,30].
All four geomagnetic stations record the changes of the geomagnetic component at the time. Figure 2 shows the long and continuous signals recorded by the four geomagnetic stations. In addition, the data record length for these stations is continuous and sufficiently long, so that there are enough time segments that can be superimposed in the Fourier transform of each period, and the data error can be well suppressed to obtain a stable C-response in the period from several days to hundreds of days. If the observation times of the station were longer than five years, the geomagnetic data of continuous records could be used to invert conductivity in the mantle transition zone [31].
When estimating C-responses, the long-term field should be removed from the highquality time series. This is because the GDS only includes the magnetic field generated by the magnetospheric current source on the earth's surface, while the long-term changing fields as the main component in the geomagnetic field observed by geomagnetic stations are generated inside the earth. Long-term fields can be calculated using the code provided by the International Association of Geomagnetism and Aeronomy (https://www.ngdc.noaa. gov/IAGA/vmod/igrf12) [32]. After that, the horizontal component of the geomagnetic field also needs to be converted from the geographic coordinate system to the geomagnetic (dipole) coordinate system. Finally, the obtained time series is processed with Bounded Influence Remote Reference Processing BIRRP software to estimate the frequency spectrum of H r and H θ corresponding to the estimated C-response according to Equation (1) and the coherence square coh 2 of each frequency, from which the data error can be estimated through the equation as ( where coh 2 is the squared coherence which can be used to assess the quality of the estimated C-responses. 1/β is the confidence level and is often set as 0.9; Z is the number of sections for a certain frequency. BIRRP is data processing software for geomagnetic sounding and the magnetotelluric response function proposed by Chave and Thomson [33]. The C-responses and squared coherence coefficients of the four stations finally obtained are shown in Figure 3. It can be seen that the C-responses at the four stations possess the same characteristics. In the short periods, the C-responses from different stations almost overlap each other, signifying that our C-responses are reliable. The errors indicated by the error bars in the short periods are rather low; in long periods, the error bars become wide, but are in the acceptable range. All these imply that our C-responses are overall of good quality. Appl. Sci. 2021, 11, x FOR PEER REVIEW 5 of 21 processed with Bounded Influence Remote Reference Processing BIRRP software to estimate the frequency spectrum of Hr and Hθ corresponding to the estimated C-response according to Equation (1) and the coherence square coh 2 of each frequency, from which the data error can be estimated through the equation as where coh 2 is the squared coherence which can be used to assess the quality of the estimated C-responses. 1/β is the confidence level and is often set as 0.9; Z is the number of sections for a certain frequency. BIRRP is data processing software for geomagnetic sounding and the magnetotelluric response function proposed by Chave and Thomson [33]. The C-responses and squared coherence coefficients of the four stations finally obtained are shown in Figure 3. It can be seen that the C-responses at the four stations possess the same characteristics. In the short periods, the C-responses from different stations almost overlap each other, signifying that our C-responses are reliable. The errors indicated by the error bars in the short periods are rather low; in long periods, the error bars become wide, but are in the acceptable range. All these imply that our C-responses are overall of good quality.  [33]. The top panel shows the square consistency; middle and bottom panels are real and imaginary parts of the C-responses.  [33]. The top panel shows the square consistency; middle and bottom panels are real and imaginary parts of the C-responses.

SA Inversion
One-dimensional inversion of GDS can be achieved by minimizing the following objective function: where λ is the regularization factor; N ω is the number of frequencies; C obs (ω) is the C-response computed from the observed data; and C mod (ω) is the model theoretical Cresponse. f m is the model constraint describing the smoothness of the model. σ is the conductivity model vector and is the parameter that needs to be inverted. In this case, the thicknesses of the thin spheres that make up the Earth are fixed and do not participate in the inversion. The process of minimizing the objective functional is similar to the process of cooling a high-energy solid material to reach a steady state of minimal energy [34]. Thus, the simulated annealing nonlinear optimization algorithm can be adapted for the one-dimensional GDS inversion. In this case, the objective function defined by Equation (3) is taken as the energy function in the annealing process. By reducing the temperature in an appropriate manner, the state corresponding to each molecule when the value of the objective function reaches the minimuml can be accepted as the solution of an inversion problem. The improved simulated annealing (ISA) algorithm proposed in this paper is based on traditional simulated annealing combined with statistical analysis to improve the robustness of inversion. The main steps in the flowchart include:

2.
Set l = 1 and perform a Markov chain loop L times.

2.1
Generate a model σ i l randomly; calculate its corresponding objective function value and RMS; 2.2 When RMS is less than RMS max , save the model as a member of the solution set for statistical analysis in step 5; 2. 3 When the objective function valueÊ corresponding to the model σ i l meets the Metropolis criterion, the model is used as the current optimal solution. The loop continues if the valueÊ corresponding to the model does not meet the Metropolis criterion; 2.4 Update l: l = l + 1;

4.
When i > M, the annealing cycle ends; otherwise, return to step 2.

5.
Statistical calculation of the optimal solution: Perform statistical analysis on all models with RMS less than RMS max to determine the final solution.
In the above iteration process, T 0 is the initial temperature; i is the current number of iterations; M is the upper limit of the number of iterations; l is the current number of cycles in the chain, which is controlled by the chain length L of the Markov chain. When the value of L increases, the number of global searches at the same temperature increases, and as a result, there is a greater probability of obtaining the global optimal solution [35]. Considering the time cost, after our experimental verifications, L = 10 can achieve the search purpose. σ max and σ min are the conductivity value boundaries, and they can be estimated based on known petrophysical models or previous experience [36]. σ 0 is the initial value of the model. SA inversion does not require a fixed initial model, so the initial value is a random number within the valid range.
RMS max represents the threshold of RMS, while RMS here is defined by An appropriate choice of RMS max is essential for inversion [25]. The selection of a reasonable RMS max needs to consider factors such as the data quality of the C-response and the number of iterations. The basic principle of RMS max setting is to require the threshold to filter out a large enough sample size in accordance with the log-normal distribution during the inversion process, otherwise the RMS max needs to be adjusted during the inversion process.
It is worth noting that in the above steps, the optimal solution obtained by the traditional SA algorithm is the final iteration result in steps (1)-(4). When step (5) is added, the inversion result is no longer a simple choice based on the result of the last iteration, but is determined by the statistical expected value of all the solutions that meet the constraint conditions during the SA inversion. Therefore, this model can better reflect the characteristics of random inversion while retaining the minimal energy state.

Random Model Generator
The random conductivity parameters conform to the temperature-related Cauchy distribution, and the specific formula is as follows: In above equations, y i is the updated factor of the i-th (i = 1, 2, . . . , N) parameter, y i ∈ [−1, 1]; T k is the temperature corresponding to the k-th iteration. u i is a uniformly distributed random variable between 0 and 1, and sgn represents a sign function. The selected and inverted results are constrained by the range of model parameters. As a result, if the model parameter interval was selected unreasonably, for example, to be in a narrower range, the inversion results would then be concentrated at the interval boundary ( Figure 4). At this time, the upper and lower limits of the parameter interval should be appropriately expanded [37,38]. An appropriate choice of RMS max is essential for inversion [25]. The selection of a reasonable RMS max needs to consider factors such as the data quality of the C-response and the number of iterations. The basic principle of RMS max setting is to require the threshold to filter out a large enough sample size in accordance with the log-normal distribution during the inversion process, otherwise the RMS max needs to be adjusted during the inversion process.
It is worth noting that in the above steps, the optimal solution obtained by the traditional SA algorithm is the final iteration result in steps (1)-(4). When step (5) is added, the inversion result is no longer a simple choice based on the result of the last iteration, but is determined by the statistical expected value of all the solutions that meet the constraint conditions during the SA inversion. Therefore, this model can better reflect the characteristics of random inversion while retaining the minimal energy state.

Random Model Generator
The random conductivity parameters conform to the temperature-related Cauchy distribution, and the specific formula is as follows: In above equations, yi is the updated factor of the i-th (i = 1, 2, ..., N) parameter, yi ∈ [−1, 1]; Tk is the temperature corresponding to the k-th iteration. ui is a uniformly distributed random variable between 0 and 1, and sgn represents a sign function. The selected and inverted results are constrained by the range of model parameters. As a result, if the model parameter interval was selected unreasonably, for example, to be in a narrower range, the inversion results would then be concentrated at the interval boundary ( Figure  4). At this time, the upper and lower limits of the parameter interval should be appropriately expanded [37,38].

Metropolis Criterion
Whether the currently randomly modified model is accepted or not needs to be judged according to Metropolis guidelines. The criterion needs to calculate the acceptance probability P of the new state by

Metropolis Criterion
Whether the currently randomly modified model is accepted or not needs to be judged according to Metropolis guidelines. The criterion needs to calculate the acceptance probability P of the new state by where ∆Ê =Ê(σ k+1 ) −Ê(σ k ) (8) represents the difference in the objective function (energy) between the modified model and the current optimal model. If ∆Ê < 0, the current solution is better than the historical optimal solution, P = 100% is agreed, and the model is fully accepted as the optimal solution in the iteration process for the next iteration. Otherwise, it means that the current solution is inferior to the historical optimal solution, and the probability P is calculated according to Formula (7) and is compared with it with a randomly generated random number uniformly distributed between 0 and 1. When the random number is less than P, the model is accepted. When the random number is greater than P, the model is not accepted and the cycle continues. The Metropolis criterion makes it possible for the solution to escape the local minimum.

Decreasing Temperature
A reasonable annealing temperature reduction mechanism can efficiently control the convergence of the iterative process to obtain an approximate optimal solution. The formula for controlling the decrease of the temperature T at the k-th iteration during the inversion process is Here, c is a constant. T 0 is the initial temperature. A high initial temperature will cause the solution to converge slowly at the beginning of the iteration, so it is not necessary to choose a high initial temperature. According to Sharma, in the DC bathymetry inversion experiment, the value of c was set to 1, the value of T 0 was set to 0.01, and the value of m was set to 2.5, which can ensure the efficiency of the inversion, and they do not depend on the model parameters during the inversion process [39]. Through experiments, we have verified that this set of parameters is also more applicable to geomagnetic depth sounding.

Statistical Analysis of the Solution Distribution
Since the conductivity parameter cannot be a negative value, all the effective model parameters in the iterative process could be described by a lognormal distribution, and lnσ obeys a normal distribution with expectation µ and variance ε 2 . For a certain layer, the final conductivity value of the inversion is equal to the expected value E(σ) of the effective conductivity parameter set of the layer, and the estimation formula is where ε is the standard deviation that can measure the uncertainty of the result. When ε is smaller, the estimated value is more reliable. The data contained within one standard deviation of the expected value account for 68% of the total data. Therefore, in this article we take this range as the reliability interval. In order to visualize the reliability of whether the samples of a conductivity parameter obey a lognormal distribution, the P-P lognormal probability diagram is a useful tool. This diagram can determine whether the sample data conform to the lognormal distribution. The P-P diagram is drawn based on the cumulative probability of the variable corresponding to the cumulative probability of the natural log-normal theoretical distribution [39]. The abscissa represents the cumulative probability density of the conductivity sample; the ordinate is the theoretical cumulative probability density, which is calculated by the probability density function When the used data obey the log-normal distribution, the cumulative probability of the sample is approximately equal to the theoretical cumulative probability, and in the P-P diagram, it will present as a 45 • straight line ( Figure 5). It can be seen that the sample distribution of the inverted conductivity σ 4 basically obeys the log-normal distribution. In the inset, the P-P diagram of the conductivity σ 4 of the fourth layer reflects this characteristic. the probability density function When the used data obey the log-normal distribution, the cumulative probability of the sample is approximately equal to the theoretical cumulative probability, and in the P-P diagram, it will present as a 45° straight line ( Figure 5). It can be seen that the sample distribution of the inverted conductivity σ4 basically obeys the log-normal distribution. In the inset, the P-P diagram of the conductivity σ4 of the fourth layer reflects this characteristic. Figure 5. The log-normal distribution histogram of the σ4 samples in the SA inversion. The inset is the P-P graph of the lognormal probability distribution of the σ4 samples. The abscissa of the subgraph is the measured cumulative probability density, and the ordinate is the expected cumulative probability density.

Results
We apply the ISA method to the one-dimensional inversion of measured geomagnetic depth sounding data from four stations around the Tarim area [40]. The above-mentioned stations are located inland at low latitudes. Therefore, their C-responses should not be affected by aurora and ocean effects [41] and thus can be directly inverted [32]. During the inversion, according to the law that the conductivity of the deep mantle jumps with the phase transition of the main minerals [42], nine layers above the core are divided [20]. The possible range for the conductivity of each layer is estimated based on the results of the Semenov inversion [22]. The control parameters of simulated annealing inversion are as follows: M = 3000; L = 10; T0 = 0.01. The control parameter RMS max is given in Figure 6. The variation range of the conductivity parameter is indicated by the shading in Figure 7. The final ISA inversion results in a smooth C-response curve with the same variation trend as the original data. At the same time, the RMS value is given in the Figure.

Results
We apply the ISA method to the one-dimensional inversion of measured geomagnetic depth sounding data from four stations around the Tarim area [40]. The above-mentioned stations are located inland at low latitudes. Therefore, their C-responses should not be affected by aurora and ocean effects [41] and thus can be directly inverted [32]. During the inversion, according to the law that the conductivity of the deep mantle jumps with the phase transition of the main minerals [42], nine layers above the core are divided [20]. The possible range for the conductivity of each layer is estimated based on the results of the Semenov inversion [22]. The control parameters of simulated annealing inversion are as follows: M = 3000; L = 10; T 0 = 0.01. The control parameter RMS max is given in Figure 6. The variation range of the conductivity parameter is indicated by the shading in Figure 7. The final ISA inversion results in a smooth C-response curve with the same variation trend as the original data. At the same time, the RMS value is given in the Figure. Figure 7 shows the final conductivity model beneath Tarim area. Above 410 km, the conductivity values below the four stations are similar. Within the sensitive depth range of geomagnetic depth sounding, there is no obvious conductivity jump under the HOQ and WEQ stations on the 410 km discontinuity, which is consistent with the conductivity characteristics of the continental transition zone obtained by Yoshino [43]. On the same discontinuity surface, the conductivity values below WUS and WMQ stations have greatly increased. We speculate that the transition zone is affected by the mantle plume beneath Tarim. From 520~660 km, the conductivity below the WMQ station in the lower layer of the transition zone is much higher than at the other three stations. At 660~900 km, compared with the lower layer of the transition zone, the conductivity of the top of the lower mantle below the four stations is significantly increased and the value is at the same order of magnitude. The noise test of the theoretical data shows that when the data contain noise, the true value of the conductivity is included in the reliable interval of the conductivity obtained by the inversion. Therefore, we believe that the inversion result of this paper is reasonable. ci. 2021, 11, x FOR PEER REVIEW 10 of 21  Figure 7 shows the final conductivity model beneath Tarim area. Above 410 km, the conductivity values below the four stations are similar. Within the sensitive depth range of geomagnetic depth sounding, there is no obvious conductivity jump under the HOQ and WEQ stations on the 410 km discontinuity, which is consistent with the conductivity characteristics of the continental transition zone obtained by Yoshino [43]. On the same discontinuity surface, the conductivity values below WUS and WMQ stations have greatly increased. We speculate that the transition zone is affected by the mantle plume beneath Tarim. From 520~660 km, the conductivity below the WMQ station in the lower layer of the transition zone is much higher than at the other three stations. At 660~900 km, compared with the lower layer of the transition zone, the conductivity of the top of the lower mantle below the four stations is significantly increased and the value is at the same order of magnitude. The noise test of the theoretical data shows that when the data contain noise, the true value of the conductivity is included in the reliable interval of the conductivity obtained by the inversion. Therefore, we believe that the inversion result of this paper is reasonable.

Model and Synthetic Data
In order to verify the feasibility and correctness of the application of simulated annealing in the one-dimensional inversion of geomagnetic depth sounding, the synthetic model was tested and the results were discussed. The theoretical model is based on Kelbert's one-dimensional model in the inset of Figure 8 [20], and its theoretical C-response (Figure 8) is inverted as observation data. The model change range during the inversion process is determined according to the maximum change interval of 400 independent one-dimensional global average models obtained by Püthe using the Monte Carlo sampling method [44].

Model and Synthetic Data
In order to verify the feasibility and correctness of the application of simulated annealing in the one-dimensional inversion of geomagnetic depth sounding, the synthetic model was tested and the results were discussed. The theoretical model is based on Kelbert's one-dimensional model in the inset of Figure 8 [20], and its theoretical C-response ( Figure 8) is inverted as observation data. The model change range during the inversion process is determined according to the maximum change interval of 400 independent onedimensional global average models obtained by Püthe using the Monte Carlo sampling method [44].

Experimental Results
We first perform an inversion test on the theoretical data without noise, shown in Figure 8. In the simulated annealing inversion process, M is set to 10,000, L is 1, and other parameters are the same as those of Sharma [38]. Figure 9 shows the poor fit and the change of model conductivity during the inversion process and the conductivity values of each layer at the end of the iteration. It can be seen that the RMS has a decreasing trend during the inversion process; since the theoretical C-response value does not contain noise and is accurate enough, according to the RMS calculation formula, its value can be reduced to a value close to zero. When the number of iterations exceeds 3000, the changes in the electrical conductivity of each layer are concentrated near the theoretical true value; when the iteration is greater than 5000, the electrical conductivity of each layer obtained by the inversion has very little difference from the true value. However, there are still cases where the poor model is accepted and the RMS is increased when using Metropolis criteria to judge the results.
Noise is widespread in actual geomagnetic data [19]. It is generally believed that the C-response affected by excessive noise cannot be used for inversion studies. Therefore, we added an acceptable 8% Gaussian noise to the theoretical C-response in Figure 8 for the inversion method test. Under the same initial conditions as the noise-free data inversion test, the inversion was performed three times using the SA algorithm. It can be seen from Figure 10a that at the end of the inversion, the C-response fitting effect of the three inversions is basically the same. However, it can be seen in Figure 10b that due to the influence of noise, the final solution (three blue dotted lines) obtained by the traditional SA method inversion is quite different. For example, the three results of σ 3 are −0.787, −0.952, −0.911 S/m. The three results are quite different from the theoretical conductivity of −0.699 S/m. It is concluded that the inversion results are affected by the difference in the sensitivity of the geomagnetic depth sounding method itself at different depths. It shows that even if the inversion is affected by noise and other factors, and even if the global inversion method is adopted, when the number of iterations is not enough, it is not easy to determine the most reasonable model. Figure 8. In the simulated annealing inversion process, M is set to 10,000, L is 1, and other parameters are the same as those of Sharma [38]. Figure 9 shows the poor fit and the change of model conductivity during the inversion process and the conductivity values of each layer at the end of the iteration. It can be seen that the RMS has a decreasing trend during the inversion process; since the theoretical C-response value does not contain noise and is accurate enough, according to the RMS calculation formula, its value can be reduced to a value close to zero. When the number of iterations exceeds 3000, the changes in the electrical conductivity of each layer are concentrated near the theoretical true value; when the iteration is greater than 5000, the electrical conductivity of each layer obtained by the inversion has very little difference from the true value. However, there are still cases where the poor model is accepted and the RMS is increased when using Metropolis criteria to judge the results. Noise is widespread in actual geomagnetic data [19]. It is generally believed that the C-response affected by excessive noise cannot be used for inversion studies. Therefore, we added an acceptable 8% Gaussian noise to the theoretical C-response in Figure 8 for the inversion method test. Under the same initial conditions as the noise-free data inversion test, the inversion was performed three times using the SA algorithm. It can be seen from Figure 10a that at the end of the inversion, the C-response fitting effect of the three inversions is basically the same. However, it can be seen in Figure 10b Figure 10 shows that when the traditional SA algorithm is used for inversion, only the last iteration value is taken as the model result, which obviously ignores the distribution characteristics of the conductivity parameters during the iteration process. Therefore,  Figure 10 shows that when the traditional SA algorithm is used for inversion, only the last iteration value is taken as the model result, which obviously ignores the distribution characteristics of the conductivity parameters during the iteration process. Therefore, if combined with statistical analysis methods, the estimated value of model conductivity may be calculated more accurately. After that, we use the method of analyzing the distribution characteristics of the solution in Section 2.3 to perform three inversions of the noise-added C-response, and the resulting model after statistical analysis is shown by the (three) red lines in Figure 10b. The initial conditions of the inversion are the same as those in Section 3. The red dotted line represents the upper and lower limits of the estimated value, which can be used as a reliable interval. Figure 10a,b shows that although the results obtained by the traditional SA method are quite different, the results of ISA are not affected by randomness, and the statistically estimated solutions are basically the same, reflecting the direction of solution convergence. The red lines in the figure basically overlap. Therefore, ISA inversion has stronger stability and accuracy, and it has weaker multi-solution. Although the data sets retained during the three inversions are different, the three final models obtained by ISA are almost identical and highly similar to the theoretical model. For example, the three results of σ 3 are −0.722, −0.736 and −0.727 S/m, which is close to the theoretical conductivity of −0.699 S/m. Although the estimated value is not exactly the same as the theoretical value, the true value of the conductivity of each layer is within the reliable interval. It should be pointed out that after experiments, the results of more ISA inversions still remain highly similar. Therefore, in actual inversion, 2~3 ISA inversions can achieve a reliable estimation of conductivity.

Reliability of the Model from Noisy Data
During the ISA inversion, the final inversion model is controlled by the choice of RMS max which is dependent on the noise level of observed data. Because our data contain noises, the reliability of our inversion model should be evaluated. In order to do so, we use synthetic data to illustrate the issue: 5%, 10%, and 12% different Gaussian noises are added to the theoretical C-response in Figure 8. These noises levels cover all possible noises in the studying data. All inversion parameters except for RMS max are the same as those in Section 3. Noise increasing demands a larger RMS max ; for the three different noise levels, after numerical tests, according to [26], the following are set: RMS 5% max = 0.08, RMS 10% max = 0.125, RMS 12% max = 0.135. Figure 11 shows the data fitting of different noises. The fitting difference obtained by the final inversion is 0.075, 0.119, 0.121 in response to the increase in the noise level. However, it can be seen from Figure 10 that the C-responses obtained by ISA inversion within the acceptable noise levels for the theoretical data are fitted very well.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 15 of 21 noises, the reliability of our inversion model should be evaluated. In order to do so, we use synthetic data to illustrate the issue: 5%, 10%, and 12% different Gaussian noises are added to the theoretical C-response in Figure 8. These noises levels cover all possible noises in the studying data. All inversion parameters except for RMS max are the same as those in Section 3. Noise increasing demands a larger RMS max ; for the three different noise levels, after numerical tests, according to [26], the following are set: RMS5% max = 0.08, RMS10% max = 0.125, RMS12% max = 0.135. Figure 11 shows the data fitting of different noises. The fitting difference obtained by the final inversion is 0.075, 0.119, 0.121 in response to the increase in the noise level. However, it can be seen from Figure 10 that the C-responses obtained by ISA inversion within the acceptable noise levels for the theoretical data are fitted very well. The inversion models are illustrated in Figure 12. Figure 12a displays the P-P diagram obtained by inversion of the conductivity of the four layers under three noise conditions. By observing the P-P diagram, we know that the effective model parameters obtained by the inversion basically conform to the lognormal distribution. Moreover, under the current noise level, the inversion result is basically not affected by noise. Figure 12b shows the final conductivity model and the reliability interval. It can be seen that although it is affected by different noises, the true values are all within the reliable range. In addition, The inversion models are illustrated in Figure 12. Figure 12a displays the P-P diagram obtained by inversion of the conductivity of the four layers under three noise conditions. By observing the P-P diagram, we know that the effective model parameters obtained by the inversion basically conform to the lognormal distribution. Moreover, under the current noise level, the inversion result is basically not affected by noise. Figure 12b shows the final conductivity model and the reliability interval. It can be seen that although it is affected by different noises, the true values are all within the reliable range. In addition, the results show that the inverted conductivity model basically coincides at 410 km below and 900 km below. On the one hand, these intervals may be in the weakly sensitive area of geomagnetic depth sounding, on the other hand, it may also reflect the ability of the method in this paper to suppress the influence of noise. In the vicinity of the depth of the transition zone, comparing the theoretical models, the inversion results better reflect the theoretical model parameters. It is worth pointing out that the noise level affects the inversion results. There is a tendency that the greater the noise, the worse the conductivity of the inversion and the larger the range of the reliability interval. Hence, when there is less than 12% noise in the data, the true model can be obtained by using the inversion method proposed by this paper. When the data noise increases again, the inverted conductivity profile will deviate greatly from the true profile, and the true profile may not be within the reliable interval. At this point, the inversion model in Figure 7 is reliable even when data contain acceptable noises.

Evidence of the Mantle Plume Beneath Tarim
At present, a few people have noticed the low velocity anomaly beneath the Tarim Basin. Kosarev pointed out that this may indicate the existence of a mantle plume. Onedimensional conductivity inversion results show that there are obvious high conductivity anomalies in the upper mantle transition zone and the upper lower mantle in Tarim area. Compared with the results of previous seismic profiles (Figure 13), the author believes that the high conductivity anomaly beneath the Tarim Basin is reliable, and the high conductivity and low wave velocity indicate the existence of a temperature anomaly.
The increase in electrical conductivity in the mantle is usually explained by an increase in the water content and/or an increase in temperature. The minerals in the mantle transition zone have a high water content, which can reach 2.0 wt% [45]. Since the significant increase in the water content in the transition zone is usually related to the dehydration of the detained plate, and there is no clear plate detention in the Tarim transition zone, it is necessary to explain the high conductivity anomaly in the mantle transition zone through the increase in temperature. In order to discuss the temperature characteristics of the transition zone in the Tarim Basin, this paper uses Yoshino et al.'s high temperature and high pressure conductivity model for the petrophysical calculation of the transition zone [43]. The conductivity profiles of 150 K, 300 K, and 650 K above the global average temperature of the mantle transition zone are obtained; σ is the volume conductivity, C w is the water content, α is the geometric factor, R is the Boltzmann constant, R = 1.38 × 10 −23 J/K, T tr is the thermodynamic temperature; σ OH , σ OP are the pre exponential factors, H is the activation enthalpy; the subscripts H and P are the electron jump conduction and proton conduction, respectively. The specific model parameters are shown in Table 2.

Evidence of the Mantle Plume Beneath Tarim
At present, a few people have noticed the low velocity anomaly beneath the Tarim Basin. Kosarev pointed out that this may indicate the existence of a mantle plume. Onedimensional conductivity inversion results show that there are obvious high conductivity anomalies in the upper mantle transition zone and the upper lower mantle in Tarim area. Compared with the results of previous seismic profiles (Figure 13), the author believes that the high conductivity anomaly beneath the Tarim Basin is reliable, and the high conductivity and low wave velocity indicate the existence of a temperature anomaly. Figure 13. Vertical cross-sections of P wave velocity perturbations along the profiles shown in Figure 1 [2]. The surface along the profile is shown at the top of the cross section. White dots show the earthquakes that occurred within a 50-km radius from each profile. The velocity perturbation scale is shown at the bottom. The dashed black lines denote the 410 km and 660 km discontinuities. The red dots are a projection along latitude.
The increase in electrical conductivity in the mantle is usually explained by an increase in the water content and/or an increase in temperature. The minerals in the mantle transition zone have a high water content, which can reach 2.0 wt% [45]. Since the significant increase in the water content in the transition zone is usually related to the dehydration of the detained plate, and there is no clear plate detention in the Tarim transition zone, it is necessary to explain the high conductivity anomaly in the mantle transition zone through the increase in temperature. In order to discuss the temperature characteristics of the transition zone in the Tarim Basin, this paper uses Yoshino et al.'s high temperature and high pressure conductivity model for the petrophysical calculation of the transition zone [43]. The conductivity profiles of 150 K, 300 K, and 650 K above the global average temperature of the mantle transition zone are obtained; σ is the volume conductivity, Cw is the water content, α is the geometric factor, ℜ is the Boltzmann constant, ℜ = 1.38 × 10 −23 J/K, Ttr is the thermodynamic temperature; σOH, σOP are the pre exponential factors, H is the activation enthalpy; the subscripts H and P are the electron jump conduction and proton conduction, respectively. The specific model parameters are shown in Table 2.  The water content in the lower mantle is extremely low, not exceeding 20 ppm [46], and there are currently very few experiments on the conductivity mechanism of the waterbearing rocks there. Therefore, when explaining the causes of high conductivity anomalies in the lower mantle beneath the Tarim area, only the temperature factor is considered. In this paper, we use Formula (13) from Sinmyo et al. [47] for petrophysical calculations for the upper layer of the lower mantle: where σ is the lower conductivity, σ 0 is the pre-referential factor, σ 0 = 74 S/m; E a and ∆V are the activation energy and activation volume, respectively, where E a = 0.7 eV, ∆V = −0.55 ± 0.01 cm 3 /mol. The values of temperature T and pressure P at different depths in the lower mantle can be found in [48]. According to Formula (12), the conductivity of the lower mantle is calculated at normal temperature, +150 K, +300 K, and +650 K, as shown in Figure 14. At the 660 km interface at the top of the lower mantle, the global average temperature is about 1925 K [43], and the electrical conductivity under the four stations around the Tarim area corresponds to the electrical conductivity after the temperature rises by 150-450 K. This shows that the temperature at the top of the lower mantle in the Tarim area has increased by at least 150 K.
bearing rocks there. Therefore, when explaining the causes of high conductivity anomalies in the lower mantle beneath the Tarim area, only the temperature factor is considered. In this paper, we use Formula (13) from Sinmyo et al. [47] for petrophysical calculations for the upper layer of the lower mantle: where σ is the lower conductivity, σ0 is the pre-referential factor, σ0 = 74 S/m; Ea and ΔV are the activation energy and activation volume, respectively, where Ea = 0.7 eV, ΔV = −0.55 ± 0.01 cm 3 /mol. The values of temperature T and pressure P at different depths in the lower mantle can be found in [48]. According to Formula (12), the conductivity of the lower mantle is calculated at normal temperature, +150 K, +300 K, and +650 K, as shown in Figure 14. At the 660 km interface at the top of the lower mantle, the global average temperature is about 1925 K [43], and the electrical conductivity under the four stations around the Tarim area corresponds to the electrical conductivity after the temperature rises by 150-450 K. This shows that the temperature at the top of the lower mantle in the Tarim area has increased by at least 150 K. Figure 14. Conductivity at 660-900 km. The global average conductivity profile is from [11]. Four orange dashed lines represent the conductivity corresponding to increasing temperatures.
It is shown in Figure 1 that the geomagnetic stations we use are distributed near the seismic profile. Figure 13 is the vertical section of the primary wave velocity disturbance corresponding to the profile. The seismic profile is from Huang's article. It can be seen that there is an obvious low velocity anomaly at 410 km below the Tarim Basin; there is a high conductivity anomaly at 410-520 km in the one-dimensional conductivity profile of WUS station, which almost coincides with the seismic profile; and there is a good corresponding relationship between the two. We speculate that the temperature of WUS station is higher than the global average of 300 K. The WEQ/HOQ station is far from the high Figure 14. Conductivity at 660-900 km. The global average conductivity profile is from [11]. Four orange dashed lines represent the conductivity corresponding to increasing temperatures.
It is shown in Figure 1 that the geomagnetic stations we use are distributed near the seismic profile. Figure 13 is the vertical section of the primary wave velocity disturbance corresponding to the profile. The seismic profile is from Huang's article. It can be seen that there is an obvious low velocity anomaly at 410 km below the Tarim Basin; there is a high conductivity anomaly at 410-520 km in the one-dimensional conductivity profile of WUS station, which almost coincides with the seismic profile; and there is a good corresponding relationship between the two. We speculate that the temperature of WUS station is higher than the global average of 300 K. The WEQ/HOQ station is far from the high temperature anomaly body accumulated at the interface of 410 km, which makes the temperature increase by only 150 K. The location of the 660 km interface penetrated by the rising hot matter may be close to the WMQ station, which is the reason why the conductivity and temperature of WMQ are higher than those at the other three stations in the range of 520-660 km. The four stations have high conductivity anomalies at the top of the lower mantle. Combined with the banded low velocity anomaly at the top of the lower mantle in Figure 13, it is speculated that due to the influence of the greater resistance of the 660 km interface, the thermal material will accumulate at the top of the lower mantle and will have a larger area, and even the thermal material will appear below the HOQ and WEQ stations.
According to the high conductivity of the upper part of the lower mantle obtained by our one-dimensional geomagnetic depth inversion, combined with the low velocity anomaly obtained by previous seismic methods, it is inferred that there is a high temperature upwelling beneath the Tarim Basin, which is hindered by the 660 km interface and accumulates at the top of the lower mantle, and upwelling occurs through the interface near WMQ and then accumulates at the 410 km interface. Moreover, the accumulation of hot material at 410 km can also be seen under Hainan volcano [3]. Therefore, the inversion results can be used as electrical evidence to speculate that there may be mantle plumes beneath the Tarim Basin.

Conclusions
In this paper, a simulated annealing algorithm combined with statistical analysis technology is used to carry out one-dimensional inversion of four stations in Tarim area. The inversion test results of theoretical data show that the simulated annealing inversion method combined with statistical analysis has stable inversion results, a strong ability to reduce noise, and the ability to restore the conductivity profile of the Tarim Basin. Based on the inversion of the measured data of four stations around Tarim, the conductivity model under the station is obtained. The results show that the overall conductivity in Tarim area is higher than the global average, especially at the top of the lower mantle and the top of the mantle transition zone. Combined with the seismic results, we believe that it is in a high temperature state in the depth range of high conductivity below Tarim. It is speculated that there may be a mantle plume beneath the Tarim Basin.