Global High-Resolution Magnetic Field Inversion using Spherical Harmonic Representation of Tesseroids as Individual Sources

: In this study, we present a novel approach combining the advantages of tesseroids in representing geophysical structures though their voxel-like discretization features with a spherical harmonic representation of the magnetic field. Modelling of the Earth lithospheric magnetic field is challenging since part of the spectra is hidden by the core field and the forward modeled field of a lithospheric magnetization is always biased by the spectral range used. In our approach, a spherical harmonic representation of the magnetic field of spherical prisms (tesseroids) is used for high-resolution magnetic inversion of lithospheric field models. The use of filtered spherical harmonic models of the magnetic field of each tesseroid ensures that the resulting field matches the spectral range of the input data. For the inversion, we use the projected gradient method. The projected gradient method easily allows us to assign an initial guess (i.e., a-priori assumption) for the inversion and avoids negative values of susceptibilities. The latter is providing more plausible models since induced magnetization is assumed to be dominant over the continents and, for the oceans, a remanence model can be subtracted. We show an application of the technique to a synthetic dataset and a satellite-derived lithospheric field model where the model geometry is based on seismic information. We also demonstrate a proof-of-concept for high-resolution tile-wise inversion for the Bangui anomaly in Africa.


Introduction
The Earth's magnetic field contains a signal from various sources including the core, ionosphere, magnetosphere, and the Earth's magnetic lithosphere [1]. The lithospheric magnetic field reflects the tectonic setting and geological provinces of the Earth, as shown by the typical stripe anomalies over the oceans reflecting seafloor spreading [2]. Often the lithospheric magnetic anomalies are interpreted from aeromagnetic surveys, where the contribution of the inducing core field is routinely reduced during processing. Aeromagnetic surveys commonly have limitations with respect to long wavelengths due to the survey extension or the methods by which individual surveys are stitched together.
Global models of the lithospheric magnetic field, on the other hand, can provide a global heterogeneous coverage and are usually based on a spherical harmonic representation [3]. For these models, the reliable spectral range is defined by the part not dominated by the main field and within the measuring bandwidth of the data sources. The main field associated with the core dominates the spherical harmonic coefficients up to 15°. Therefore, lithospheric models like LCS-1 [4] only define a reliable lithospheric field spectral content for more than 15°.
LCS-1 provides the coefficients to degree and order 185 based on measurements from the CHAMP and Swarm satellite missions. Reference [5] showed that a combination of the gradient observations with the vector information improves the lithospheric field determination, and the reliable range can be increased to degree n = 154 in comparison to only vector measurements, which can reach degree n = 138 only. However, for previous models like the various generations of MF (e.g., [3]), it is stated that satellite-data provide the full signal to degree and order 85, while some signal loss occurs at higher degrees. LCS-1 is using data from the multi-satellite mission Swarm and shifts this range to 133° and provides even higher coefficients. This improvement is significant since the reliable spectral range of the satellite lithospheric field is currently one of the basic limitations concerning the robustness of interpreting the lithospheric field model in terms of magnetization and magnetic crustal thickness on continental or global scales.
There are experimental data sets derived from Swarm (and older CHAMP) data, which are nonregularized based on the inversion approach from References [6] and [7]. These already shift resolution up to 200° and should be used for further future analysis.
A number of studies previously presented models of the lithospheric magnetization within the spherical harmonic degree range corresponding to the satellite-derived lithospheric field. First magnetic crustal models were presented by References [8] and [9], based on seismic crustal thickness estimated and susceptibility estimations from geology. These models took in account only an induced magnetization and were not coinciding with Magsat anomaly maps [10]. The Standard Earth Magnetization Model [11] included both induced and remanent magnetizations, using the 3SMAC model [12] to set different tectonic areas (i.e., magnetic properties) and a remanent magnetization in oceans from Reference [13]. This model was updated several times by adding further crustal models and information about heat-flow [14]. References [15,16] introduced a model where ranges of susceptibilities and magnetizations were assigned for each geological region of the world, according to the geographical information system (GIS). The model was adjusted via iterative inversion to fit satellite-derived (mostly from CHAMP satellite mission) lithospheric field model MF7 [3]. The result was the vertically integrated susceptibility (VIS) model and vertically integrated remanent magnetization (VIM) for the oceanic crust. Later, the remanent model was modified by making plate reconstruction corrections and improving oceanic remanence estimates [17]. However, the VIS model is still considered the reference for any later attempts.
Another limitation in inferring information for the magnetic sources comes from the nonuniqueness of the magnetization solution due to the existence of magnetization structures that do not produce any visible magnetic field, which are known as annihilators [18][19][20]. That means that any inversion requires a priori information or a reference model to result in plausible susceptibilities.
All these inversion and modelling approaches used fitting and/or inversion based on point dipole sources. That allows only a simplified representation of the geometry of the crust, especially when seismic and other types of information are considered. However, modelling of aeromagnetic surveys requires a better representation of the source bodies and, therefore, often Cartesian or Fourier-based approaches are used [21]. A representation with magnetic tesseroids [22] provides a versatile discretization that can be used for modelling both satellite and airborne data. The forward modelling code of Reference [22] allows the calculation of the magnetic field in spherical coordinates using spherical prisms, which are also called tesseroids (magnetic sources).
However, magnetic inversion of satellite-measured data within a particular wavelength window, corresponding to the reliable lithospheric field, is only viable if the signal of the individual sources is also filtered to the same spectral range. Therefore, the calculated field (across the entire spherical surface) of each tesseroid can be converted into a spherical harmonic representation. This enables the inversion for the reliable wavelength range while keeping the geometrically confined and discretized lithospheric magnetization model. This set-up allows a reliable calculation of the magnetic field at any altitude above the Earth's surface.
In the following, we first outline the inversion method that is designed to have an adequate source discretization and inverts for a specific spherical harmonic range only. The validity of the approach is demonstrated with a synthetic example and next applied to the satellite-model LCS-1 in order to provide a global lithospheric susceptibility model, which, as we demonstrate at the end, can serve as a reference or background model for further studies on local scales.

Tesseroid Forward Field
The software package magnetic tesseroids [22] is used for the forward calculation. Magnetic tesseroids calculates the vector magnetic field of a given tesseroid model on a specified grid. The magnetic contributions of all tesseroids in the model are summed up in the resulting data grid of the magnetic field vector components. In case of only induced magnetization, the vector magnetic field of each tesseroid is directly proportional to its magnetic susceptibility [22]. Where , , and are the rotation matrices.
P is the tesseroid's gravity tensor in the coordinate system of the computation point P and C is the magnetizing field in the coordinate system of the tesseroid's center C.
Therefore, the magnetic field vector in each point of the two-dimensional spherical grid can be expressed by the equation below.
where = 1 … , = 1 … are the indices of the computational point within this spherical grid with longitudes and latitudes.

Spherical Harmonic Model
In this study, we utilize the software package SHTools [23] to convert the forward calculated spherical grids ( P ( , ) ) into a spherical harmonic representation. For the inversion, we set calculated grids ( P ( , ) ) to be equidistant with the number of points = 2 × . To convert the grid into a spherical harmonic model, SHTools utilizes the algorithm of efficient computation of the spherical harmonic expansion of functions defined on a two-dimensional sphere [24]. The relative errors using the method on = 2 × sampled grids are similar to grids with datapoints at Gauss-Legendre quadrature spacing (see Reference [24] for more details).

Projected Gradient Inversion
Generally, geophysical inversion is ill-posed. Due to the presence of noise in the observed data and rounding errors in forward calculation, the existence, uniqueness, and stability of the inversion are hard to guarantee. For unique solution, constraints or a specific inversion method must be applied. However, in classical methods like least-squares inversion or Tikhonov regularization, the resulting susceptibilities are not constrained. Hence, the inversion can result in negative susceptibilities, which are unreasonable if occurring for large areas over the continents. In case of a constrained inversion over a small-scale anomaly, a strong misfit between the observed field and forward calculated field of the inversion result would indicate possible remanent magnetization. However, large-scale remanent sources are uncommon for the continental crust and for the scale of our analysis.
Therefore, we adopt the projected gradient method [25], which allows us to exclude negative susceptibility values and to set a priori constraints.
The projected gradient method is a way to solve a constrained optimization problem, or, more specifically, a non-negative least square linear inversion problem. The optimization problem is illustrated in the formula below.
where is the sensitivity matrix, which represents the magnetic field spherical harmonic coefficients of each tesseroid, is the magnetic susceptibility, and represents the spherical harmonic coefficients of the observed magnetic field model.
The cost function ( ) is shown below.
where = , = 2 , and the constrained optimization problem becomes a constrained quadratic programming (QP) problem.
The Gradient Descent type method is a standard way to solve the quadratic optimization problem. To ensure a non-negative solution to the optimization problem, the projection operator is used to force the iteration back to the constraint range. Specifically, the projected gradient method iteratively uses the gradient projection to update the current solution to +1 by the following rule. where is the step size and ( ) is the search direction which is defined by the equation below.
To find a reasonable step size in each iteration, there are sub-iterations to either increase or decrease the initial guess 0 = 1 by 0.1 until satisfies the sufficient decrease condition.
However, the computational complexity per iteration is relatively large to find a step size satisfying the condition. To reduce the cost of checking this condition, Lin (2007) uses a Taylor expansion of the cost function and simplifies the sufficient decrease condition as shown in the equation below.
We use the Projected Gradient Method with the simplified form of the sufficient decrease condition for our analysis, in the way it is presented here.

Data
For both the synthetic and real data application, we use a similar set of data to confine the geometry and parameters space. For the synthetic inversion, a global lithospheric field model is forward calculated from a crustal model with predefined susceptibility. The same geometry is later used for the inversion of a measured lithospheric magnetic field model.

Model Geometry
For all calculations, we assume the entire crust to be magnetic and set the geometry of the model, which defines the size of the tesseroids. The top of the model is the top of the crystalline basement as provided from the model CRUST 1.0 [26], originally with a 1° × 1° resolution (Figure 1a), but, for the global application, it is resampled to a 2° × 2° resolution. The base of the model is the Moho boundary. The Moho depth model is taken from Reference [27] (Figure 1b). The model is based on a similar catalogue as the Moho depth in CRUST 1.0, but, based on a geostatistical kriging approach, it provides uncertainties, which can be potentially exploited in further modifications. This set-up defines a one-layered crust, which is presented with tesseroids with variable thickness.

Magnetization Direction and Susceptibility
Induced magnetization in the crust is a response of the mainly ferrimagnetic rock properties to the inducing main field. We use International Geomagnetic Reference Field model, year 2014 [28], to define the direction and strength of the inducing field in the center-point of each tesseroid.
For the initial susceptibility distribution in the model, we use the VIS model [16]. For each tesseroid, the susceptibility is calculated from their model by dividing the VIS model by the height of the tesseroid. The VIS value used is the mid-tesseroid value (Figure 2a).
To explore the effect of a-priori information (initial guess) on the inversion, we calculate a simplified susceptibility distribution with averaged values for oceanic and continental crust ( Figure  2b). We use the map of the age of the oceanic crust [2], separate the oceanic and continental crust, and calculate the mean average susceptibility from the VIS model in Figure 2. This results in an average value for the continental crust of 0.0154 SI and for the oceanic crust of 0.0555 SI.
While we neglect remanence in the synthetic tests, we also consider remanent magnetization as a priori information in the application to real data. For this, we use the remanent magnetization vectors from Reference [17] and calculate the corresponding field using 1-km thin tesseroids to imitate original dipole methods of Reference [17]. The model of Reference [17] relies on many initial assumptions, such as shape of the sources in the oceanic crust, which may not coincide with the actual sources in the real crust. However, it is currently the best available model to subtract the remanent field of the oceans.

LCS-1
For the global inversion of real data, we choose the lithospheric magnetic field model LCS-1, which is the spherical harmonic model of the lithospheric field derived from the CHAMP and Swarm satellite observations (Olsen et al., 2017). Figure 3a shows the field for spherical harmonic degrees 16-89 at 400-km height, which represents satellite orbit height. The lower cut-off is in line with the original model, while the higher cut-off is chosen to be in line with our model parametrization.  We use the remanent magnetization model of Masterton et al. (2013) [19] and calculate the corresponding magnetic field for the spherical harmonic degrees . For this, we assign the remanent vector model to a spherical shell of 1-m thickness made of tesseroids 0.25° in longitudinal and latitudinal width. From this, we calculate the remanent field (Bz component) at the altitude of 400 km (Figure 4). Truncated spherical harmonic model of the calculated remanent field is shown in Figure 4b and the residual after subtracting this field from LCS-1 in Figure 3b.

Synthetic Observed Field
For the analysis, we create a synthetic input data set for testing our inverse approach. For this, we use magnetic tesseroids [22] to calculate the magnetic field of a crustal model at the satellite altitude of 400 km. The grid spacing is chosen to be 1°, which approximately corresponds to a spherical harmonic degree of n = 89, due to RAM limitations. The calculated grid is converted into a spherical harmonic model with a maximal degree N = 89 using SHTools [23] (see Figure 5a). To further imitate a satellite-derived lithospheric magnetic field model, we truncate all degrees below and including degree n = 15. The filtered synthetic observed Bz component is shown in Figure 5b.

Synthetic Test
The first test explores how the inversion of a truncated field recovers the initial true model used to create the synthetic field. The loss of the long wavelength part between 1° and 15° implies that the inversion might tend to result in solutions that may not match susceptibility distribution in the true model.
The projected gradient method requires an initial solution that is used as a first guess before iterations start. To explore how the truncation of the field affects the inversion results, we introduce two a priori models: x 0 with the susceptibility distribution from the true model (see Figure 2a), and x 0 with averaged susceptibility values for oceanic and continental crustal areas (Figure 2b).

A-Priori Model from a True Susceptibility Model
This first test helps identify the error introduced by grid processing and creation of spherical harmonic models of each tesseroid field within the algorithm. For the inversion, using a truncated field with a true model 0 as an initial assumption, we expect that manipulations with grids introduce a rounding error and a complete fit cannot be achieved. After the first iteration, the projected gradient norm reaches its minimum (approx. 75,696 × 10 −4 ) and remains the same for all further iterations.
One can use the iteration value of 75,696 × 10 −4 for the linear norm to set up a threshold for future inversions, since this value represents the error introduced by re-gridding and represents the conversion to spherical harmonics. RMS of the difference between the true model H 0 and the inversion result H N=1 after the first iteration is about 0.0 SI, i.e., the non-zero linear norm value after the initial iteration does not affect the result and the way floating point numbers from the result are stored in the output text files (max 18 digits after the point).

Initial Guess with Averaged Susceptibility for Oceanic and Continental Crust
Next, we explore the effect of the a priori model. Earth lithosphere magnetization is not well known. Hence, the effect of the a priori assumptions might seriously affect the results. For this, we compare the results of the inversion (Figure 6b) with the initial guess x A 0 (Figure 2b) with the true model x H 0 (Figure 2a). Figure 7 and Figure 8 indicate that the fit for the field is very good for degrees n > 15 and the RMS of the difference between the observed and calculated field is very low (0.00933 nT) (Figure 6b). However, the long-wavelength part of the spectra is very different. The resulting susceptibility distribution A N is similar in short-wavelength patterns to the true model x H 0 (Figure 2a), but has an overall decreased amplitude, especially for spatially large areas (Figure 7a). This can be explained by the differences in the long-wavelength range. The inversion provides only susceptibilities in the limited spectral range. To explain degrees 16-85, the amplitude of long-wavelength anomalies (and, therefore, susceptibility of spatially wide sources) can be lower than in the true model. That is also seen in the spectral representation in Figure 8.

Inversion of LCS-1
For the inversion with real data, we again choose VIS-based model H 0 . The geometry of the model is the same as in the synthetic example and we only replace the synthetic data set with the observed lithospheric magnetic field. For the inversion, we invert LCS-1 directly after removing the remanent field as an additional apriori parameter.
The inversion result is shown in Figure 9a, and the comparison of the result's spectrum against the original LCS-1 is shown in Figure 9b. RMS of the difference between the observed truncated field and forward calculated truncated field of the result LCS N=10,000 is 0.03938 nT.

Performance of Inversion
A synthetic inversion case with averaged susceptibilities for oceanic and continental crust as an initial guess ( A 0 ) showed that powers of spherical harmonic degrees for N = 1000 iterations are, overall, lower than for N = 10,000 iterations (see Figure 8b). A minimal value for the linear norm of projected gradient after N = 10,000 iteration is approximately 2157.838 × 10 −4 (iteration 9999), which is two magnitudes higher than the threshold estimated in the very first test.
A higher number of iterations appears to increase the number or artifacts in the inversion result (see Figure 6 for example), which can be attributed to the existence of short-wavelength anomalies (with the wavelength close to approximately 1°) in the input that cannot be caused by 2° × 2° tesseroids. Thus, the issue is a regularization problem, which is not discussed here.

Use of a-priori Model
Our results confirm that a good a-priori model for long-wavelength anomalies is necessary for accurate results. Currently, the existing dataset that can be used as an initial guess for the inversion of real data is a Hemant VIS-based model H 0 .
Since synthetic inversion for the initial guess, which is different from the true model, showed deviation in long-wavelength below degree n = 16, one cannot expect that long-wavelength sources were recovered correctly in the inversion result (Figure 9). A similar conclusion was made in Reference [29]. This is especially evident for the inversion results where the initial model in Reference [16] does neither represent the detailed surface geology nor a realistic susceptibility model due to the absence of such information (e.g., parts of Africa and Antarctica). The inversion does change the initial guess and leads to a further refinement of the sources, but, apparently, the amplitude stays as low as possible to fit the short-wavelength.
The application to the Bangui anomaly ( Figure 10 and Figure 11) shows how the inversion changes the initial susceptibility model by adjusting the local values. Clearly, the choice of the a priori information largely influences the inversion results. Therefore, defining the a priori model with uncertainties could lead to further improvement by adding a stochastic element to the inversion. However, such a model does not exist on a global scale.
Therefore, an alternative would be to use actual rock susceptibility measurements and their statistical properties. Databases of such a measurement exist for Norway [30], Fennoscandia [31], and Canada [32] and can be used to provide an initial guess of a mean susceptibility and its standard deviation. A refined inversion could be done for those areas to explore the influence by a priori information with more detail.

Layering of the Model
Another shortcoming of our current models is that the vertical structure of the crust is not well expressed by a single layer. Measurements from outcrops may not adequately represent the entire crust an since the susceptibility of rocks beneath the surface might be different. Therefore, further developments could be adding more information on the crustal structure from seismic information or models like CRUST 1.0. For example, Reference [16] have based their model on a classification of the ratio of magnetization in the upper and lower crust. The lower crust is, in general, assumed to be less magnetic than the upper crust. However, the most prominent satellite lithospheric anomalies can be expected to have a significant contribution from the lower crust or even the upper mantle [33]. Therefore, one of the next steps should be to add the possibility of layering the crust and to modify the inversion for this. While this is already possible with the current set-up, one has to consider depthweighting or regularization or the use of aeromagnetic data in the inversion.

Regularization
The other solution to help the inversion resolve long-wavelength sources more realistically is to make a stochastic inversion, which would control lower degree power spectra at each iteration. Such an approach would still require an initial guess for power spectra of lower degrees (such as Reference [34]), derived from the field of an a-priori initial susceptibility guess/model. Control over power spectra would result in the appearance of long-wavelength anomalies in the result, which would correspond to highlighted spherical harmonic degrees.
Such measures would impose implicit regularization on the inversion. Explicit regularization in the inversion, such as a way to guarantee smoothness of the result, may be detrimental to fitting a shorter wavelength. Figure 6 shows that the N = 10,000 solution, which is less smooth than the N = 1000 solution, fits the short-wavelength spectrum better. We do not discuss explicit regularization in the inversion further because we expect to use the tile-wise approach later to refine the result.

Application to Airborne Data -Proof of Concept
While there are many potential ways to enhance the analysis, the current model can already be used as a reference in modelling regional data. For areas on a regional or local scale, the current model can be used for the surrounding areas to avoid edge effects, but might also serve as a first guess for the model domain.
A tile-wise approach can be applied, when we select a moving window-tile and calculate the field from the areas surrounding this tile using susceptibilities from the global inversion. Then, we can subtract this far-field effect from the observed field and repeat the inversion for this tile with higher resolution. Usage of tesseroids as sources instead of dipoles allows scalability of the method, and this approach can be used for the inversion of airborne data. In addition, inversion results can be tested against airborne data, which can be helpful to refine inversion results and models in the future.
As a proof of concept, we demonstrate the results of the tile-wise inversion with grid and tesseroid model resolution of 0.5 degrees for the Bangui Anomaly area. We have chosen this area because it features the most prominent anomaly at satellite height. Although remanence was proposed as a primary source for this anomaly (see Reference [35]), that is not necessarily true and depends on the geophysical data used as a constraint. For the demonstration of the concept, the nature of the anomaly is secondary.
We used data without a remanent field from Figure 4b, and removed the forward calculated surrounding field as an input. x LCS N=10,000 is used as an initial guess for tile inversion. The forward calculated field of the inversion result at the airborne altitude (5 km) is shown in Figure 10 against the original LCS-1 (up to 179°). One can see that the field on the airborne altitude was recovered through the inversion of the data on the satellite altitude, even though there are artifacts related to discretization (Figure 10b).   Figure 11 shows the process of refinement of the results for Bangui area. Figure 11a demonstrates the initial assumption H 0 , while Figure 11b shows the global inversion result with LCS-1 data (from Section 4.2.). Lastly, Figure 11c shows the tile inversion for this area with the resulting Figure 11b as an initial assumption.

Conclusions
We developed an approach for global high-resolution inversion of satellite magnetic data using spherical harmonic models of tesseroids. Our method utilizes algorithms from the SHTools [23] package to convert forward calculated grids of tesseroids into spherical harmonic models and Projected Gradient [25] to perform the inversion.
Our approach is compared to classical spherical harmonic approaches that are very versatile and easily allows us to couple different geophysical, voxel-based approaches. For example, it can be easily linked to the LitMod3D software [36], which is a popular software for integrated geophysical and petrological modelling of the lithosphere.
Using the developed approach, we investigated caveats of the global magnetic inversion. One major finding was that selection of a-priori information, i.e., susceptibility distribution in the crust, significantly affect the long-wavelength part of the inversion results for a forward calculated magnetic field.
There are two possible ways to ensure more accurate inversion results. The first one is to enforce a particular power for lower degrees via stochastic inversion in order to fit some suggested lowdegree power spectra. The other option is to use measured susceptibilities across the globe to set up a more accurate a-priori model. However, the global model presented in this paper can already serve as a background model for more detailed local modelling studies.