Geological Mapping of Grib Kimberlite Pipe Using 3D Magnetization Vector Inversion From Airborne Data

1. Introduction

Grib kimberlite pipe is located in the Archangelsk diamondiferous province (north of European Russia). It was discovered in 1996 and is now one of major diamond producing mines in Russia [1]. It has estimated reserves of 98.5 million carats of diamonds and an annual production capacity of 3.62 million carats. Its surface area is measured more than 16 hectares and is traced by drilling operations to the depth of 1024 m. Since its discovery, the pipe has been a subject to multiple remote sensing and geophysical surveys, photogrammetry as well as extensive drilling operations, which allowed constructing an accurate geological model. All remote sensing and geophysical data sets collected since 1996 were in some way subject to cultural and/or industrial disturbance [2]. The study of the pipe itself is further complicated by the overlaying 60-80 m thick Cretaceous sediments, which make its mapping impossible using satellite remote sensing and photogrammetry techniques and challenging even with usage of airborne magnetometry.

Generally, kimberlitic pipes, which are the main sources of diamonds consist of three facies: root (deeper) facies, diatreme facies, and crater (upper) facies [3]. In many cases, the pipes are eroded to some extent, and usually the diatreme facies, composed by hypabyssal kimberlite with typical magmatic structure and composition, emerges at the surface level. However, in some areas, kimberlite with lesser degree of erosion may expose bowl-shaped crater zone at the surface, formed by a variety of eruptive and resedimentation processes. The kimberlitic material of the crater facies is often diluted with sedimentary material, making it difficult to identify by direct observations, remote sensing and geophysical methods.

1.1. Geology

Grib pipe is located in central part of Zimniy bereg (Winter shore) diamondiferous province within Chernoozersky kimberlite field. It was formed during D₂-C₁ magmatism, cutting through Rifean and Vendian sedimentary strata and is overlain by terrigenous-carbonate formation of middle Cretaceous age, which, in turn, is overlain by loose quaternary sediments (Figure 1).

Rifean sysem is represented by interbedded silts and clays with thin lenses of variegated sandstones. Its estimated thickness is 50 – 70 m.

The upper Vendian system is composed of four suites: Ust-Pinezhskaya (V₂up); Erginskaya (V₂er) Melskaya (V₂ml) and Zolotnickaya (V₂zl). V₂up is characterized by predominance of silts over sandstones with presence of three tuffaceous horizons. Estimated thickness of V₂up pack is 560 m. V₂er has thickness of 64 to 87 m while V₂ml has thickness of 165 to 180 m. Both are composed of interbedded silts and sandstones. Perhaps the thickest pack of sediments is V₂zl reaches up to 400 m in thickness and is composed of red-coloured fine-grained sandstones with interbeds of silts.

Middle Cretaceous and Quaternary rocks are directly overlaying the kimberlite pipe with Cretaceous system reaching 52 to 82 m and comprised of uniform quartz sandstones slightly mixed with gravel material at the base of the pack.

Quaternary system is composed of sands, sandy loams and clays and reaches from 0.4 to about 20 m in thickness. There are unstratified dolomites and dolomitic limestones present directly below quaternary pack varying in thickness from 10 to 30 m.

The kimberlitic rocks of Grib pipe can be conditionally stratified into crater and hypabyssal facies. Crater facies is the shallowest part of the kimberlite with depths interval from 52 to 209 m (averaging 100 m in thickness). It is composed of various rocks including sandstones mixed with igneous material, tuffaceous sandstones, tuffs and sedimentary breccias.

Hypabissal facies can be followed to depths exceeding 1000 m below surface and is composed of two types of material. The upper part of hypabyssal facies is composed of kimberlitic tuff breccias, the other type is pure kimberlite. Thin carboniferous veins are widely present. The contact of kimberlite with host rocks are sharp with abrupt and clear transitions [1].

1.2. Magnetic properties

This section is mainly based on research of paleomagnetic properties of the rocks composing the Grib pipe, as well as the host rocks present within the survey area carried out by VNIGRI Russian Research Institute of Oil and Mineral exploration [4].

The magnetic properties of rock types from area under study have a very wide range. (Table A1) contains scalar characteristics of J_n (magnitude of natural remanent magnetization measured in A/m); K_m (magnetic susceptibility in SI units); Q (Koeniksberger ratio); Inclination (I_o); Declination (D_o); statistical parameter describing likelihood that the values lie outside of 95% confidence range (α95), sample masses and direction of natural remanent magnetization of kimberlites and kimberlitic tuff breccias.

The histograms of J_n and K_m distribution are shown in Figure 2. Table 1 contains mean values and ranges of natural remanent magnetization and magnetic susceptibilities of kimberlite rock types from Grib pipe, Table 2 contains magnetic properties of host rocks from area under study.

The relationships between scalar values of Jn with Km are shown in Figure 3. They form elongated groups with visible linear trend, which reflects concentration of magnetic minerals in the rock types. The latter may be attributed to potential preservation of primary paleomagnetic information. Some data from kimberlites and tuff breccias overlap.

The solid lines in Figure 3 represent averaged values forming linear trends.

Evaluation of induced magnetization (Ji) can be done without accounting for demagnetization coefficient due to body shape, as the values of magnetic susceptibility are low (Km <<0.1 SI units). Induced magnetization can be derived for each rock type using the mean values of Km from Table A1 and parameters of current magnetic field (HT) in the area of study (D = 17.9°; I = 77.3°; HT = 43.6 А/м at 65.5°N/41.4°E), Ji = Km*HT. The values of Ji are also provided in Table A1 on last lines of every rock type section and comprises 22% of the corresponding Jn values for kimberlites and 87% for kimberlitic tuff breccias and 24% for all other rock types.

Considering wide spread of sample collection locations, Ji does not have significant effect on the area in general, however there is a narrow localization of significantly increased Ji values, following trend of samples 44, 45, 49, 50, 51. 66, 67 and 68 (Figure 4), which generally corresponds to central occurrence of type 2 (consolidated kimberlite) hypabyssal facies of Grib pipe.

Summarizing the materials presented in this section (more comprehensively studied in [4]) the following generalized statements can be made:

• Sampled rocks have the following average values: J_n = 2.32Е-01 А/m, K_m = 1.29E-03 SI units. а. J_i=43.6∙1.29E-03=5.61E-02 А/m.
• The highest average values recorded in kimberlites are as following: J_n = 2.80E-01 А/m, K_m= 1.41E-03 SI units, J_i = 43.6∙1.41E-03 = 6.13E-02 А/m.
• Consolidated kimberlites form anomalous group of samples with high anisotropy of magnetic susceptibility (P = 1.58). These samples are characterized by high values of J_n, reaching 5 A/m. On the map of Grib pipe these samples fall within a compact area elongated in NW-SE direction.
• The distribution of J_n/K_m is best described by linear relationships, which, in turn, reflects concentration of magnetic minerals in sampled rock types.
• The predominant magnetic minerals responsible for J_n are iron hydroxides (magnetite, titanomagnetite and hematite).
• The directions of full natural remanent magnetization have a wide range. Both direct and reversed polarity is observed.
• The average direction of induced magnetization of sampled rock types has been evaluated using anisotropy of magnetic susceptibility (D=27°, I=77°), which is consistent with direction of current magnetic field in the area under study, J_i = 43.6∙1.29E-03 = 5.61E-02 А/m

2. Materials and Methods

We made a comprehensive study of an airborne data set, which was collected in 1985, some 11 years prior to kimberlite pipe discovery. The airborne magnetic data set was aquired using a fixed-wing Antonov-2 aircraft with average terrain clearance of 120 m. There was no obvious response from the pipe in the 1995 airborne data, therefore the data were first subject to trend removal, then subject to two types of 3D geophysical inversions:

• Magnetization vector inversion [5]
• QDIK inversion for susceptibility and remanent magnetization [6]

Both resulting data sets were compared to known geology as well as to each other in order to establish correlation with kimberlite facies uncovered by drilling, as well as to suggest an improved routine for advanced interpretation of airborne magnetic data for kimberlite exploration in cases when remanent magnetization is present.

The total magnetic intensity data collected over Grib kimberlite pipe in 1985 are shown over bedrock geology in Figure 5а. The flight lines are trending EW and are rendered in black colour. As it can be seen from the figure, there is no obvious reflection of the pipe in the RMF data, therefore the data were subject to trend removal Figure 5b ,which resulted in some improvement in mapping the kimberlite, however with only fractional reflection of the target.

Methods of 3D magnetic data inversions described in this section include full magnetization vector inversions: MVI [5] and QDIK [6], which both allow to recover magnitude and direction of full magnetization vector J_t = J_i + J_n.

3. Results

3.1. MVI inversion

The abbreviation stands for Magnetization Vector Inversion [5]. This is a voxel (3D) inversion carried out on a regular rectangular mesh, which allows to recover both remanent and inductive magnetization by means of solving inverse problem for strength and direction of full natural magnetization from magnetic intensity data and without prior knowledge of the direction or strength of remanent magnetization. In majority of magnetic inversion algorithms it is assumed that the magnetic field is entirely a result of magnetic induction. In many field data sets however there are strong and obvious effects of remanent magnetization, including cases, when total magnetic intensity yields negative values.

Magnetic Induction B(r_j) is a vector field, which can be described in a form of equation (1) in magnetostatics [7]

$$B\left(r_j\right)=\nabla {\int }_{V}M\left(r\right)·\nabla \frac{1}{\left|r-r_j\right|}d{r}^3$$

(1)

where r_j is the observation point; V is the volume containing magnetization vector M(r). In case of MVI inversion we invert B(r_j) for M(r) using Tikhonov minimum gradient regularizer (2):

$$\begin{gathered} Min\varphi \left(m\right)={\varphi }_D\left(m\right)+\lambda {\varphi }_M\left(m\right) \\ {\varphi }_D\left(m\right)=\sum _j^M{\left|\frac{G_{j}m-B_j}{e_j}\right|}^2 \\ {\varphi }_M\left(m\right)=\sum _Y^3{\left|w_Y{\partial }_{Y}m\right|}^2+{\left|w_{0}m\right|}^2 \\ \lambda :{\varphi }_D\left(m\right)=X_T^2 \end{gathered}$$

(2)

where 𝜙(m) is the total objective function, composed of data objective 𝜙D(m) and model objective function 𝜙M(m). The two parts of the total objective function are scaled using the trade-off parameter λ, which is chosen based on convergence of the routine in chi-squared sense χ². The data objective function includes the forward modelling operator G_jm, which is at each iteration step is compared with observed data B in order to fit within the assigned error e_j. The model objective function contains gradient of the model Δ_γm and the altitude of the model with weighting terms w_γ and w₀. In this routine additional constraints such as upper and lower bounds can be placed on m, which can be derived from a-priori information.

The TMI data collected over Grib pipe was subject to MVI inversion. The data set was optimized to accommodate one sample per mesh cell. A regular mesh measuring 208 x 210 x 153 cells in its core area was used for the MVI inversion with smallest cell measuring 15 x 15 x 10 m in X Y and Z directions accordingly. Padding of 665 m (6 cells starting from 32m m with expansion coefficient of 1.5) was allowed in every direction. Convergence to target misfit was achieved in 16 iterations).

As a result, several voxel models were calculated including amplitude of magnetization vector. Corrected magnetic susceptibility values were calculated from magnetization amplitudes by using a multiplication factor of 43.6 (derived from the petrophysical studies). The susceptibility voxel model derived from MVI inversion is shown in Figure 6. Magnetization vector voxel model is shown in Figure 7. The data fit of MVI inversion is shown in Figure 8.

3.2. QDIK inversion

This is an example of so-called “cascade” inversion used on potential field data. The methodology described in this section is applied to magnetic data and describes recovery of both inductive and remanent components of magnetization vector from magnetic survey data by means of inversion using adaptive meshes [5,8,9]. This is a multiparametric inversion capable of recovering magnetic susceptibilities alongside (K) with Koeniksberger ratio (Q) and remanent Inclination (I) and Declination (D).

The possibility of global minimization of the residual field at a consecutive allocation of local objects, which provides a monotonous reduction of dispersion of the residual field is shown. For the inversion of the field in the local area, the algorithm of robust regression by the method of elastic network with the attraction of prior information on the physical properties of geological objects is used. Results of the inversion of the magnetic field allow defining magnetic susceptibility, and vector of residual magnetization of local objects, and density for the gravitational field.

The described algorithm uses different types of regularization of inverse problem, which allows accurate recovery of total magnetization vector unique for each mesh element. Regularization algorithms vary in boundary conditions and weighting functions incorporated in model objective function [9].

In order to separate inductive component of magnetization vector from remanent it is necessary to know the distribution of magnetic susceptibility over the area of modelling, therefore it is important to use a-priori information, which may be defined in form of other inversion results, range of predicted magnetic properties and possible shapes of the magnetized targets.

The following section describes the algorithm of the QDIK-inversion in greater details

Input data

• Magnetic field data, corner coordinates, geometric and physical parameters of inversion models:

  -

  Anomalous magnetic $∆\mathbf{T}$);

  -

  Digital elevation model and sensor height;

  -

  Normal magnetic field parameters ( $T_0$– magnitude of normal magnetic field vector; $D_0$– declination; $I_0$– inclination);

  -

  A-priori petrophysical information as ranges of possible physical property values:

  ◦

  Magnetic susceptibility (${\mathsf{\chi }}_{min}$, ${\mathsf{\chi }}_{max}$ );

  ◦

  Remanent magnetization vector parameters and their discretization: ($∆\mathrm{*}$):

  ▪

  Koeniksberger ratio ($Q_{min}$, ${\mathrm{Q}}_{max}$, $∆\mathrm{Q}$);

  ▪

  Declination (${D_r}_{min}$, ${D_r}_{max}$ , $∆\mathrm{D}$);

  ▪

  Inclination (${I_r}_{min}$, ${I_r}_{max}$, $∆I$);

  -

  Discretization of media:

  ◦

  Mesh geometry: regular discretization of the media volume under study and padding distance to compensate for out of the mesh sources of magnetic fields;

  ◦

  Dimensions of the rectangular mesh elements.

• Algorithms and parameters used to select magnetic anomaly contours and corresponding local areas of the mesh:

  -

  Median Absolute Deviation (MAD) function used to evaluate field anomaly in sliding window mode; standard deviation function and other statistical tools;

  -

  Rectangular sliding window range of sizes along corresponding coordinate axes: ${L_x}_{min}$, ${L_x}_{max}$; ${L_y}_{min}$, ${L_y}_{max}$ ;

  -

  Measure of field anomalous properties set as threshold to determine presence of local:  ${\epsilon }_{min}$;

  -

  The parameters allowing to evaluate sizes of mesh model for selected local anomaly and the respective area of the field for local inversion.

Generalized inversion workflow

А: Detailing looping of inversion models

Definition of geometric parameters and petrophysical constraints for the inversion model.

Definition of statistical functions and parameters for local target selection

QDIK-inversion of the magnetic field of j-th mesh model

${∆T}_{mod}=∆T;$${∆T}_j={∆T}_{mod}$; ${\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}_{{dT}_j}={‖{∆T}_{mod}‖}_{L_2}$B: local field anomalies inversion loop ${∆\mathit{T}}_{\mathit{j}}$

• Local anomaly contour selection subroutine $C_L$, alongside with corresponding area of fragmented mesh model $V_L$ and inverted ${∆T}_L$area;
• Formation of $p$ elements of possible remanent magnetization parameters multitude $QDI=\{Q,D_r,I_r\}$, acquired from ranged sequences based on parameter ranges: ($Q_{min}$, ${\mathrm{Q}}_{max}$, $∆\mathrm{Q}$); (${D_r}_{min}$, ${D_r}_{max}$ , $∆\mathrm{D}$); (${I_r}_{min}$, ${I_r}_{max}$ , $∆I$);

C: Vector scanning loop: $\mathit{k}=1,\dots ,\mathit{p}$QDIK-inversion of magnetic field ${∆T}_L$ based on ELNET algorithm [6,10] with definition of $(Q_k,{D_r}_k,{I_r}_k)$ parameter selection, responsible for minimization of model objective function ${{∆T}_{mod}}_L$ with data misfit ${∆T}_L$.

Here ${∆Tmod}_L={{\delta T}_{mod}}_L+b_L$, where ${{\delta T}_{mod}}_L$- is anomalous effect; $b_L$- background value for the local inversion area.

C: End loop

3.

Decision on local object selection as a result of changing norm of initial field

if ${‖{∆T}_j-{\delta T}_{mod}(Q_k,{D_r}_k,{I_r}_k)‖}_{L_2}<{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}_{{dT}_j}$; $\delta Tmod$ – anomalous effect from a local object, calculated for the entire range of field definition.

then

${∆T}_j={∆T}_j-{\delta T}_{mod}(Q_k,{D_r}_k,{I_r}_k)$;

${\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}_{{dT}_j}={‖{∆T}_j‖}_{L_2}$;

Blocking of mesh elements for selected local object;

Unblocking of previously blocked elements of local areas.

else

blocking of local area elements to prevent infinite looping.

B: End loop: absence of new selected anomalous areas or free mesh elements for local object selection.

Analysis of local magnetic object selection results and possible correction of selected object multitude by means of inversion. The total number of selected objects is $N_{obj}$.

D: Selection and correction of  ${\mathit{N}}_{\mathit{o}\mathit{b}\mathit{j}}$parameter loop

Selection of k-th object from $N_{obj}$ responsible for maximum correlation of its anomalous effect with residua field ${∆T}_j={∆T}_j-\sum _{i\ne k}^{N_{obj}}{{\delta T}_{mod}}_i$;

E: Layer extension loop for k-th object to ensure its shape (9 iterations by default)

Forming of $p_k$ elements of possible parameter value multitude for remanent magnetization ${QDI}_k=\{Q,D_r,I_r\}$, derived from ranked rows, based on corrections achieved in the D-loop for the following ranges: ($Q_a$, ${\mathrm{Q}}_b$, ${∆\mathrm{Q}}_{ab}$); (${D_r}_a$, ${D_r}_b$ , ${∆\mathrm{D}}_{ab}$); (${I_r}_a$, ${I_r}_b$ , ${∆I}_{ab}$).

F: Vector scanning loop: $\mathit{k}=1,\dots ,{\mathit{p}}_{\mathit{k}}$QDIK-inversion of ${∆T}_L$ field based on the ELNET algorithm with determination of $(Q_k,{D_r}_k,{I_r}_k)$ prameters, responsible for minimization of model objective function ${{∆T}_{mod}}_L$ and data misfit in local area ${∆T}_L$.

Here ${{∆T}_{mod}}_L={{\delta T}_{mod}}_L+b_L$, where ${{\delta T}_{mod}}_L$- anomalous effect; $b_L$- background value of the local inversion area.

F: End loop

E: End loop with premature termination in case of data misfit increase.

A decision is made to update local target parameters, based on perturbations to the observed fields

if ${‖{∆T}_j-{\delta T}_{mod}(Q_k,{D_r}_k,{I_r}_k)‖}_{L_2}<{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}_{{dT}_j}$; ${\delta T}_{mod}$ – anomalous effect of the local object, calculated for the entire area of magnetic field definition.

Then${∆T}_j={∆T}_j-{\delta T}_{mod}(Q_k,{D_r}_k,{I_r}_k)$;

${\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}_{{dT}_j}={‖{∆T}_j‖}_{L_2}$;

Blocking of mesh elements of the accepted local object.

else

exclusion of k-th object from the list of selected objects and unblocking of corresponding mesh elements of the model.

D: End loop: in absence of new selected anomalous areas or unblocked free cells for new local target selection.

А: End loop: decision of the user to end loop manually and/or to advance to a new model included within the are covered by the observed data.

This algorithm was applied to the airborne magnetic data measured over Grib pipe. The starting mesh elements were set to 400 m in every direction and refined three times to 25 m. Figure 9 shows resulting model with all non-zero elements as well as data fit. A voxel model of magnetic susceptibility values was constructed from the QDIK inversion results (Figure 10).

4. Discussion

The results of two 3D inversion attempts are presented in a side by side comparison with known bedrock geology. Figure 11, Figure 12, Figure 13 and Figure 14 show plan view slices at different absolute elevations.

It becomes evident that both MVI and QDIK inversions show anomalous susceptibilities corresponding to pipe contours. The results are further examined in a cross-section. Figure 15 shows the AB geological cross-section over the Grib pipe. Figure 16 show sampled voxel models imposed over this cross-section (a): MVI inversion results and (b) QDIK inversion results.

Finally, the result of the inversions were compared with the direct core petrophysical measurements. Figure 17, Figure 18, Figure 19 and Figure 20 show comparison of susceptibilities Koeniksberger ratios and values of remanent Inclination (I_r) and Declination (D_r). Voxel models presented in these figures were sliced at absolute elevation of 6 m, which corresponds to approximately 100 m depth (the upper part of the kimberlite pipe).

5. Conclusions

In current research we have demonstrated two different approaches to recover three-dimensional distribution of magnetic properties. It is evident that the MVI inversion yields smooth, yet confident recovery of the consolidated kimberlite units at depths starting at 100 m from surface (Figure 16), which is consistent with petrophysical studies. Furthermore, MVI approach allows noninvasive corrected magnetic susceptibility recovery from airborne data, which is at the same time somewhat consistent with susceptibilities uncovered by detailed sample petrophysical studies.

The QDIK inversion is quite complicated in terms of practical application however it provides interpretable results, including corrected susceptibility, Koeniksberger ratio and remanent magnetization vector parameters, all recovered non-invasively through analysis of remotely measured (airborne) data.

Comparing the two approaches, it can be concluded that MVI algorithm is more suitable for large-scale data sets and should be adapted as standard workflow in cases when remanent magnetization effects can not be neglected, while the QDIK algorithm is preferred on target-scale inverse modelling of non-invasive (airborne) data sets.

This study provides workflow for interpretation of both large and local airborne magnetic data sets, which is seen as a powerful tool in geological mapping, kimberlite exploration and the only algorithm existing to date to increase geophysical recovery of otherwise blind targets. It can be further combined with other remote sensing data sets to provide more effective target recovery.