CC BY
11
RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 21, ES3002, doi:10.2205/2021ES000766, 2021
Numerical algorithms for structural magnetometry inverse problem solving
P. S. Martyshko1, N. V. Fedorova1, and A. L. Rublev1
Received 22 March 2021; accepted 20 April 2021; published 21 May 2021.
New methods for studying the anomalous magnetic field of the lithosphere and modeling the sources of anomalies in different layers of the Earth's crust are given in this article. The interpretation consists of three stages: separation of anomalies from the sources in different layers of the Earth's crust, data conversion to the magnetic pole and solving the inverse problem of finding the surface of sources. To isolate anomalies from different layers of the earth's crust, we used a technique based on the subsequent continuation of magnetic data up and down. Calculation of the vertical component of the anomalous magnetic field by its absolute value and conversion to the pole is based on the approximation of the anomalies of the magnetic induction modulus through a set of singular sources - magnetized rods. To solve the nonlinear inverse problem and calculate the boundaries separating layers with different magnetization, the method of local corrections was applied. The new computer technology was created using parallel computing in the multiprocessor computer systems. To demonstrate the developed methods, we chose the region of Western Urals 300 km x 230 km, where the large magnetic anomalies are located. Anomalies from various layers of the Earth's crust were identified and the models of the sources of magnetic anomalies in the granite layer and the surface of the basaltic layer were constructed. The study revealed basic-ultramafic massifs under a thick sedimentary cover, as well as extended deep belts and large elevations of the basalt layer of the Earth's crust. The results clarify the lithosphere structure of the Western Urals and can be used for geodynamic reconstructions. KEYWORDS: Magnetic anomaly; Earth's crust; modelling; parallel algorithms; inverse problem.
Citation: Martyshko, P. S., N. V. Fedorova, and A. L. Rublev (2021), Numerical algorithms for structural magnetometry inverse problem solving, Russ. J. Earth. Sci., 21, ES3002, doi:10.2205/2021ES000766.
1. Introduction lithosphere and building geomagnetic models of the
Earth's crust.
Magnetic anomalies provide important informa- Magnetic surveys using satellites have allowed
tion about the structure and the history of the the pattern of the magnetic field distribution to
Earth's lithosphere. Study and analysis of the be acquired and the centennial variations of the
anomalous field structure have practical value for field to be specified over the whole globe, with-
geological mapping, geological survey, and prospect- out "white spots". At present, a knowledge about
ing for mineral resources, and in the sense of reveal- the spatial-temporal peculiarities of the terrestrial
ing the peculiarities of the tectonic structure of the magnetic field has been significantly broadened,
—^--—— Tn t~> a i_ii , and the methods of spherical harmonic analysis
institute of Geophysics UB RAS, Ekaterinburg, , , ,, , , , ™ . , „ ,
t} have been well developed. lhis has allowed sci-Russia
entists to elaborate the international models of the
Copyright 2021 by the Geophysical Center RAS. main magnetic fields (International Geomagnetic
http://rjes.wdcb.ru/doi/2021ES000766-res.html Reference Field or briefly IGRF) for every 5-year
ES3002
1 of 9
Figure 1. Tectonic maps Eastern Europe (a) and Western Urals region (b), Russia. Uralian megazones: W - West, C - Central, T - Tagilo-Magnitogorskian. 1 - Western boundary Preuralian foredeep. 2 - Uralian deep faults and their numbers: 1 - West, 2 -Central, 3 - Main Uralian.
interval using the unified technique [IAGA Division, 1995].
In the framework of the Federal Special-Purpose Program "Development of prognostic-geophysical maps for the principal mineralgenic zones of Russia", the territory of the Russian Federation has been mapped with digital maps in the past decade. The information basis for the maps were the data of areal airborne and land-based magnetic surveys of large, middle, and small scales made. The geophysical data were processed using modern computer methods. In order to unify the survey data and unite them, observations on the reference airborne magnetic profiles were made and the Uralian mapping network was organized. The correction values of centennial variations in the geomagnetic field were calculated based on the IGRF [IAGA Division, 1995; Olsen, 2002]. Resulting from these works, a digital model of the map for the anomalous geomagnetic field ATa in the Uralian Region and its vicinity (East European Plain and West Siberia) was created [Chursin et al., 2008]. We used a set of this data to model the sources of magnetic anomalies in the layers of the Earth's crust.
We have developed new methods and applied an efficient computer technology to study the structure of the anomalous magnetic field and to build the sources of anomalies in different layers of the Earth's crust. In the study of large anomalies, it is necessary to process big data, therefore, in our new computer technology we used parallel computing on multiprocessor computer systems. The results of using this method are demonstrated for the Western Ural region, covering the Eastern part of the ancient East-European platform and the Western part of the Ural orogeny 300 km x 230 km (Figure 1).
The paper is organized as follows: we present mathematical algorithms for separating magnetic field sources in depth, for the approximating of anomalous magnetic field by array of magnetized rods and for 3D magnetic data inversion. Then all algorithms are applied to the Western Urals (Russia) to build the sources of anomalies in different layers of the Earth's crust (to extract the effect of the interfaces between layers and find its topography).
2. Methods
2.1. Method of Separation of Anomalies From Sources in the Earth's Crust
To divide the long- and short-wave components of the amplitude spectrum of anomalies, geophysi-cists utilize numerical methods of field simulation at various altitudes [Martyshko and Prutkin, 2003; Martyshko et al., 2014].
The anomalous magnetic field has an integral character and contains components from all the sources located in the upper lithosphere. In order to extract the anomaly from the sources in different layers of the Earth's crust, a technique based on subsequent upward and downward magnetic data continuation was used [Martyshko et al., 2014]. There are two problems to be solved. Firstly, we continue the observed data upward to the height H to diminish the influence of the sources in the near-surface layer (up to a depth H). Then, we continue the obtained function downward to the depth H, i.e. to the distance 2H in downward direction.
Since the problem of downward continuation is a linear ill-posed inverse problem, therefore we must use some regularization.
Finally, the field was continued upwards again to the Earth's surface. So, we receive the magnetic field of the sources that are located below the depth H. Subtracting of this field from the observed one gives the field of the sources that are located within the layer (from Earth's surface to the depth H). By repeating this procedure in different heights and depths, we could separate the fields generated by the layers between the corresponding depths.
Let [x, y, z} are the rectangular Cartesian coordinates with axis Z pointing downwards and plane XOY coinciding with the surface of the observations and ATa = U. Recalculation of the magnetic field U(x, y, z)lz=o measured on the Earth's surface area upward to the level z = —H is made by the Poisson's formula:
U(x,y, —H ) =
2k
' —OO 'J — oo
[(x - X1 )2 + (y - y')2 + H2]3/2
Fredholm equation
1 rr
2k
' —oo J — oo
2H
[(x - x')2 + (y - y')2 + 4H2]3/2
U(x', y', H)dx'dy' = U(x, y, -H), (2)
where U(x',y',H) are field values at the depth H to be determined.
Solving the integral equation of the first kind
(2) is an ill-posed problem, which requires application of the regularization methods. The operator of equation (2) is positive definite and self-adjoint; therefore we can apply the scheme of Lavrentiev [1967].
After the grid discretization of equation (2) and approximation of the integral operator by the quadrature formulas, the problem is reduced to solving the system of linear algebraic equations (SLAEs) with the symmetric matrix K. In the regularized form, the system is
(K + aI)U = U, (3)
where U is the field values on the depth H, I is the identity matrix and a is the regularization parameter.
For solving the SLAEs (3) the iteratively regularized simple iteration method was used
uk+1 = uk-
1
Am
(K + al) uk - U , (4)
where Amax is the maximal eigenvalue of matrix K + al (symmetric case), k is the number of iteration.
The stopping criterion in the iterative processes
(4) is the fulfillment of the condition
Kuk - UI
< e
U(x',y', 0)dx'dy'. (1)
Recalculation of the magnetic field U(x, y, —H) to the depth H is implemented by the first type of
\\U ||
with a sufficiently small e.
Parallel algorithms. At the stage of upward recalculation of the magnetic field during calculation of integrated operator (1), the parallel algorithm of the matrix to vector multiplying is used. Paralleling of the gradient type iteration methods
x
is based on splitting the K matrix by vertical bands into m blocks, and the solution vectors u and the vectors U from the right part of the system of algebraic equations into m segments, in such a way that n = m x L, where n is the dimension of the equation system and m is the number of processors.
At any iteration, each of m cores (processors) computes its part of the solution vector. If the K matrix is multiplied by the u vector, each one of m cores (processors) multiplies only the respective part of lines from the K matrix to the u vector. The host processor performs data transfer and also computes its part of the solution vector.
The parallel algorithms of the height transforms based on the gradient type iterative methods have been designed and numerically implemented on Uran supercomputer at the Institute of Mathematics and Mechanics of the Ural Branch of the Russian Academy of Sciences [Martyshko et al., 2014]. The Uran supercomputer is a cluster of multiprocessor systems, each computing unit of which is based on two 4-core Intel Quad-Core Xeon processors with a speed of 3 GHz and a 16/32 Gb cache. Overall, 1664 computing cores and 3584 Gb cache memory are available to the user.
We note that for solving the problem by the double-precision iterative methods on the Uran supercomputer, significant memory is required for storing matrix K. For the 300 x 300 grid, the number of the elements of matrix K takes 60.3 Gb, while the code-run system on the Uran cluster normally allocates at most 2 Gb for each separate process of the MPI program. Therefore, matrix K is preliminarily divided by the vertical bands into blocks and is stored in the memory of each processor by parts. Thus, considering the memory constraints, at least 50 Uran processors are required to solve our problem.
The computational experiments show that for large data volumes, paralleling the algorithms significantly reduces the time of the computations [Martyshko et al., 2014].
2.2. The Approximation of Anomalous Magnetic Field by Array of Magnetized Rods
The magnetic surveys are mainly carrying out on the equipment that measured the absolute values of the geomagnetic field T. So the anomalous mag-
netic field of the lithosphere is often represented by ATa. Since, ATa is not a harmonic function, for correct application of methods for magnetic inverse problems solving it is necessary to determine the vertical component Z. The calculation method for the vertical component of an anomalous magnetic field based on its absolute value was developed. Conversion is based on the approximation of magnetic induction module anomalies through a set of singular sources and the subsequent calculation for the vertical component of the field with the chosen distribution. In a paper [Mayer et al., 1985], it is proposed to use a set of rods that are uniformly magnetized along a geomagnetic field. Each rod is determined by seven parameters: six coordinates of two edge and magnetization of the source. Magnetic anomalies from such sources are expressed through algebraic dependencies and are calculated very quickly. It was shown that a small number of such sources can well approximate complicated magnetic anomalies [Mayer et al., 1985].
The approximation reduces to solving the following nonlinear programming problem: minimize the functional, which is the difference between the measured data and the field of model sources:
M
F(N,p) = ^(ATa, - f (N,p,xl,yl,zl))2, (5)
i= 1
where ATai is the field values specified at the points (xi, yi, Zi) on observation surface (may be in irregular grid); N is the number of objects involved in the procedure; P is the vector of rod parameters with dimension 7 x N; f is the function that calculates the total field value of N rods, parameters of which are determined by the vector P at the points
Model rods are small parametric of source class, it increases stability of the inverse problem solution. To minimize the functional (5), various nonlinear programming methods were used and the Polack-Ribier method was selected, which converges faster than the others [Byzov et al., 2017]. The algorithm was implemented using the parallel computing technology on the NVIDIA GPU (GeForce GTX 780).
The method allows high-precision approximation of complex magnetic anomalies. Using the solution of the direct problem for the found model rods, it is easy to calculate the values of Z for the vertical
y
(x, y, 0) X
-^
Ii (x, y, z(x,y))
H
z V h
Figure 2. Two-layer model of the Earth's crust.
direction of the magnetization vector. For anomalies, the intensity of which varies from 1000 to 3000 nT, the error of the obtained values is estimated as ±20 nT. Note the conversion to the magnetic pole can significantly reduce the computational process for solving the inverse problem to find the geometry of the volume sources.
2.3. Method to Solve the Inverse Problem
To solve the equation (6) and find the function z(x,y) we use algorithm based on the modified method of local corrections [Martyshko et al 2010]. Method of local corrections was proposed for the approximate solution of nonlinear inverse gravimetric problems [Prutkin, 1986]. It is based on the assumption that the variation of the field at a certain point is affected mostly by the variation of the part of the object boundary closest to this point.
A solution is calculated from data automatically by successive approximations, without a time-consuming trial-and-error process. A contrast value of physical properties and the position of a horizontal reference plane should be specified to obtain a unique solution.
Discretization of equation (6) leads to the following system of nonlinear equations:
c Kioio (zv) = Uio jo > (7)
Let {x, y, z} are the rectangular Cartesian coordinates with axis Z pointing downwards and plane XOY coinciding with the surface of the observations. We consider the two-layer medium in 3D space. The model consists of two layers of a constant magnetization. Surface separating the upper and lower layers at a sufficient distance from the center of area goes to the asymptote (Figure 2).
We calculated the vertical component of the magnetic induction Z of the contact surface, separated the layers with vertical magnetization, I\ and I2, at the point (x,y, 0) on the ground surface by the formula:
where c is the weighting factor of the cuba-ture formula, Ui0j0 AZ(xí0 ,yj0, 0) is the left-hand side of equation (6), Zij = z(xí,yj,), Ki0j0 = K(xí0 ,yj0 ,X{,yj ,Zij) is the integrand. The scheme of an iterative method of finding the borders of magnetized layers defined by the equation z = z(x,y) is as follows. At each step of solving (7) we try to reduce the difference between the given and the approximate values of the field at a given node solely by modifying the value of the unknown function at that node.
As a result, we received an iterative formula for finding z^1:
/+œ r+œ
/ x
-œ J—œ
(
z (x,y )
[(s - ¿)2 + (y - y')2 + 32(X,y)]3/2
H
[(x - x')2 + (y - y')2 + H2]
3/2
)
=
V )
1 + «(4)2 x (Ul3 - un)'
(8)
dxdy, (6)
where z = z(x,y) is the equation of the surface separating the upper and lower layers, the AI is the magnetization jump at the boundary layers, H is the horizontal asymptote.
where a is the regularization parameter, {z™j} are the values of the unknown function z(x,y), n is the number of iteration. The values of the regularization parameter in (8) are selected based on the results of numerical experiments. In this case a = 1E — 3. At lower values a, the process converges to the same solution, but more slowly. At large values a, the process diverges.
(
Figure 3. Magnetic anomalies of the Western Ural region (a) and separated magnetic anomalies from upper layer of the Earth's crust 0-5 km (b), granite layer 5-20 km (c) and basalt layer (d). The designations are shown in Figure 1.
3. Case Study: Sources of Magnetic Anomalies of the Western Urals
To demonstrate the developed methods, we chose the region of Western Urals 300 km x 230 km, within which large magnetic anomalies are located (Figure 1a).
In tectonic terms, the territory covers the ancient Precambrian East European platform, and the western structures of the Ural fold system. The Uralian orogen is located along the western flank of a huge (> 4000 km long) intracontinen-tal Ural-Mongolian mobile belt. The orogen developed mainly between the Late Devonian and the Late Permian [Puchkov, 2009]. The Urals are divided into several north-south striking structural zones, giving the Urals a general appearance of an approximately linear fold belt. The following structures are located within our site: the Preuralian fore deep, the West Uralian mega zone, the Central Uralian mega zone, the Tagilo-
Magnitogorskian mega zone. Structures are divided by deep faults: Western, Central and Main Uralian (Figure 1b).
The area is well studied by seismic methods; a generalization of the research results is given in the monograph [Druzhinin et al., 2015]. The thickness of the Earth's crust within the East European platform reaches 40-45 km and increases to 55 km under the Uralian orogen. According to seismic data, the Earth's crust consists of three main layers: sedimentary, granite and basalt. The thickness of the sedimentary cover on the East European platform and the Preuralian fore deep exceeds 4 km, and within the Urals, the foundation is brought to the surface of the earth. The boundary between the granite and basalt layer is on average located at a depth of about 18-20 km and can vary from 1430 km.
The intensity of the magnetic anomalies reaches 1000-1500 nT (Figure 3a). Digital data from the map of the Ural region were used to study the structure of the anomalous magnetic field [ Chursin
Figure 4. The results of approximation of AT anomalies (sources - rods and the depth of their edges are shown) (a) and reduction to the magnetic pole (b). Z map is colored, and AT contours are shown with dashed lines.
et al., 2008]. The anomalies from magnetized massifs in the upper layer of the Earth's crust to the depth of 5 km were identified using up-down transformations for the conversion height H = 5 km and were obtained by subtracting these data from the initial field (Figure 3b). For magnetized sources located in deeper layers below 5 km, up-down calculations were performed for H = 20 km. Anomalies from magnetized blocks in the granite layer were obtained as the difference between the calculated fields for H = 5 and 20 km (Figure 3c). The magnetic field calculated for H = 20 km consists of long-wave regional anomalies that correspond to the integral distribution of magnetization in the lower basalt layer of the Earth's crust (Figure 3d).
For the separated magnetic anomalies of the module induction AT, approximation was performed by singular sources - rods magnetized along the modern geomagnetic field. For the regional component of the magnetic field, the approximation results and the segment's locations, as well as their depths are shown in Figure 4a. Under the epicenter zones of the anomalies, the depths to the sources are minimal and amount to 16 km. Then, the anomalous fields were reduced to the magnetic pole and, using the direct problem, the values of the vertical component Z were calculated for the found parameters of the sources with vertical magnetization. The difference between the anomalies of the module induction AT and Z is shown in Figure 4b. Despite the fact that our region is located at high latitudes of 560 — 580 and the geomagnetic field vector has an inclination close to vertical 730,
the epicenters of the Z anomalies are shifted 1215 km north of the epicenters of the total induction anomalies (Figure 4b).
To solve the inverse problem, the magnetic field data from each layer were specified on grids of 100 x 100 points. As a result of an iterative process, the geometry of the sources in each layer was determined, while the relative deviation error of the initial and model anomalies was less than 13 percent, with the exception of the edges of the region of 10-20 km. Therefore, we reduced the region and cut our models 20 km from each edge. When choosing the source magnetization in the layers, we based on the results of previously performed calculations along seismic profiles, as well as three-dimensional models of regional anomalies [Fedorova, 2005; Fedorova and Rublev, 2019; Fe-dorova et al., 2017]. The main carriers of magnetization in the Earth's crust are the minerals of the titanomagnetite series and, above all, magnetite. As a rule, these minerals are found in significant quantities only in basic and ultramafic rocks, which can have a high magnetization of 2-6 A/m. Sedimentary rocks contain a small amount of magnetic minerals and this layer does not create noticeable magnetic anomalies. The average magnetization of the granite layer is low, not exceeding 0.3 A/m. According to geothermal estimates, in our region, at depths of more than 30-35 km, the temperature exceeds 5800 (Curie point of magnetite) and the minerals lose their ferromagnetic magnetization. Therefore, it is possible to limit the magnetic sources in a basalt layer 30 km deep. Models of
Figure 5. The results of modeling sources in upper 0-5 km layer of the Earth's crust (a), the granite layer (d) and the surface of the basalt layer (c). The designations are shown in Figure 1.
sources of the magnetic anomalies in the granite and basalt layers were built for magnetization of 3 and 4 A/m.
The modeling results in the form of a three-dimensional surface of the sources in the layers are shown in Figure 5.
In the upper layer (Figure 5a), small sources are present within the East European platform but their depth is more than 4.5 km. In the Urals, east of the Main Uralian Fault, magnetized massifs are located close to the surface of the Earth, which is in good agreement with the data of geological mapping. To the west the Tagilo-Magnitogorskian megazone are bordered by ser-
pentinitic melange. This megazone predominantly consists of Ordovician-Lower Carboniferous complexes of oceanic crust and island arc, including the Platinum-bearing Belt of layered basic-ultramafic massifs, overlain by platformal carbonate and rift-related volcanic rocks [Puchkov, 2009].
In the granite layer, magnetized sources form linear belts 100-200 km long (Figure 5b). Probably, on the East European platform, these are green belts formed at different Precambrian times. The strike of the zones varies from sub-latitudinal to submeridional, which indicates a long and complex history of the formation of the ancient platform. In the Urals, 3 belts are extended in the north-south
direction. They correspond to the deep-Uralian faults, which extend along the Ural structures for more than 3000 km.
The basalt layer has 5 large elevations in the form of domes, which rise up to 15-17 km, and the surface deflections reach 26 km. (Figure 5c). Note that the basaltic layer domes are present not only within the East European platform, but also under the Urals structures up to the Main Uralian Fault. Consequently, the upper layer of 15-20 km of the Earth's crust of the Urals within the Western and Central megazones was formed as a result of thrust from the east and lies on the ancient basaltic layer of the Precambrian platform.
We note that values of the regularization parameter in equation (3) depend on downward depth. For example, if depth is 5 km then a = 0.23, if depth is 20 km then a = 0.208.
4. Conclusion
With a help of the new computer technology, the structural features of the anomalous magnetic field for the Western Urals territory were studied. Anomalies from various layers of the Earth's crust were identified and the models of the sources of magnetic anomalies in the granite layer and the surface of the basaltic layer were constructed. The study revealed basic-ultramafic massifs under a thick sedimentary cover, as well as extended deep belts and large elevations of the basalt layer of the Earth's crust. The results clarify the lithosphere structure of the Western Urals and can be used for geodynamic reconstructions.
Acknowledgments. The Russian Science Foundation (Grant No. 20-17-00058) is acknowledged for supporting this research.
References
Byzov, D., L. Muravyev, N. Fedorova (2017), The approximation of anomalous magnetic field by array of magnetized rods, AIP Conference Proceedings, 1863, 560051, Crossref Chursin, A. V., A. M. Prut'yan, N. V. Fedorova (2008), Digital anomalous magnetic field map of the Northern, Middle, Southern Urals and adjacent territories of the East-European and West-Siberian platforms, Litosfera, 6, 63-72. (in Russian)
Druzhinin, V. S., P. S. Martyshko, et al. (2015), Tectonodynamic Model of the Earth's Crust of the Urals and Neighboring Areas, Doklady Earth Sciences, 463, No. 1, 659-662, Crossref Fedorova, N. V. (2005), Modeling of the magnetic field variation in relation to the origin of the Manchazh secular trend anomaly, Izvestiya. Physics of the Solid Earth, 41, No. 5, 355-362. Fedorova, N. V., A. L. Rublev (2019), Numerical Modeling of Magnetic Anomalies Sources in the South Urals Earth Crust, Russian Geology and Geophysics, 60, No. 11, 1310-1318, Crossref Fedorova, N. V., A. L. Rublev, et al. (2017), Magnetic anomalies and model of the magnetization in the Earth's crust of circumpolar and polar the sectors of Ural region, Geofizicheskiy Zhurnal, 39, No. 1, 111122, Crossref IAGA (1995), Division V, Working Group 8, International Geomagnetic Reference Field, 1995 Revision, J. Geomag. Geoelectr., 47, No. 12, 1257-1261, Crossref
Lavrentiev, M. M. (1967), Some Improperly Posed Problems of Mathematical Physics, 76 pp. Springer, New York. Crossref Martyshko, P. S., I. L. Prutkin (2003), Technology for separating the gravity sources by the depth, Ge-ofiz. Zh., 25, No. 3, 159-170, (in Russian)Cross-ref
Martyshko, P. S., N. V. Fedorova, et al. (2014),
Studying the Structural Features of the Lithospheric Magnetic and Gravity Fields with the Use of Parallel Algorithms, Izvestiya. Physics of the Solid Earth, 50, No. 4, 508-513, Crossref Martyshko, P. S., I. V. Ladovskii, A. G. Tsidaev (2010), Construction of regional geophysical models based on the joint interpretation of gravitaty and seismic data, Izvestiya. Physics of the Solid Earth, 46, No. 11, 931-942, Crossref Mayer, V. I., F. I. Nikonova, N. V. Fedorova (1985), A numerical optimization procedure for interpreting gravity and magnetic-anomalies, Izvestiya, Fizika Zemli, 5, 46-57. Olsen, N. (2002), A model of geomagnetic field and its secular variation for epoch 2000 estimated from Orsted data, Geophysical Journal International, 149, No. 2, 454-462, Crossref Prutkin, I. L. (1986), The solution of three-dimensional inverse gravimetric problem in the class of contact surfaces by the method of local corrections, Physics of the Solid Earth, 22, No. 1, 49-55. Puchkov, V. N. (2009), The evolution of the Uralian orogen, Geological Society, London, Special Publications, 327, No. 1, 161-195, Crossref
Corresponding author:
P. S. Martyshko, Institute of Geophysics UB RAS, 620016, Ekaterinburg, Russia. (pmart3@mail.ru)
