Seismic Density Model of the White Sea’s Crust

: Study of the deep structure of the White Sea region is relevant to active geodynamics, manifestations of kimberlite magmatism, and the prospects of oil and gas searches. The aim of this work was to model the velocity and density structure of the earth’s crust in the White Sea region. Modelling was carried out using the known data of instrumental observations and the software complex “Integro”. With the help of 2D models based on deep seismic sounding (DSS) proﬁles and a digital map of the anomalous gravity ﬁeld, the density structures of local areas of the region’s crust were reﬁned. A 3D density model was built. Within the framework of this model, the positions of the density layers were determined. The relief of the Mohorovichich (Moho or M) discontinuity reﬂects the anomalies of the gravity ﬁeld. Depression of the Moho boundary in the bottleneck of the White Sea indicates the vertical structure of the earth’s crust associated with manifestations of kimberlite magmatism.


Introduction
The subject of the current study is the White Sea basin and adjacent territories. Located at the junction of two large tectonic elements of the East European Craton, the Fennoscandian Shield and the Russian Plate, this region is constantly experiencing dynamic loads caused by the continuing uplift of the Fennoscandian Shield. Its original crustal structures formed in the Archean were partially transformed in the processes of Proterozoic rifting and subsequent tectonomagmatic activation. Studies of geodynamics, tectonics, and the evolution of the material composition of the lithosphere are relevant in the region. Its characteristic feature is the manifestation of kimberlite magmatism, and deposits of diamonds and other minerals. It is believed that the Arkhangelsk province, which ranks second in Russia in diamond mining after Yakutia, is far from exhausting its diamond potential. The geological research recently conducted in the region aimed to find hydrocarbons. The formulation and solution of theoretical and applied problems are facilitated by the study of the deep structure of the region.
Initially, geological and geophysical study of the White Sea region [1] was based on the results of aeromagnetic and gravimetric surveys and drilling of wells [2]. Using the analysis of geophysical fields and seismic data, the positions and connections of the structural elements of the Earth's crust of the region were clarified, and the idea of the White Sea Rift System of the Rifeian emergence was formed [3]. A feature of this system was the absence of an anomalous mantle under the riffs, characteristic of modern formations. Riftogenic deflections, gradually filled with sedimentary material, have traditionally been considered possible oil and gas collectors [4].
Interpretation of geological and geophysical data is usually carried out in conjunction with the construction of 2D and 3D models using the petrological characteristics of rocks, such as density and magnetization. Several geological and geophysical models of the deep structure are known for the White Sea region [15][16][17][18][19]. Many of the models derived from references [15][16][17][18][19] are characterized by detail loss, incomplete information, uneven coverage of areas, and differences in local volumes of data used. Some [17,19] are calculated by solving the direct problem, representing the geological and geophysical environment in cylindrical blocks with selected sizes and positions, and searching for the solutions which better fit the measured fields. This approach is currently being revised with modern modelling methods and tools, which allow solving direct and inverse geophysical problems without resorting to idealizing the shapes and boundaries of geological objects.
One such tools is the Integro software package developed by VNIIgeosystem for solving predictive and diagnostic problems and problems of thematic regionalization of territories [20]. The analytical apparatus of the complex automates the solution of direct and inverse problems of geophysics, allowing digital maps to be drawn, cartographic references to be carried out, and 3D data to be processed, visualized and stored. The algorithms of the complex work with regular and irregular data networks, and allow spatial analysis of vectors and surfaces to be conducted, and necessary sequences and cross-sections of objects to be built. The resulting models are distinguished by their detail; these models visualize geological structures and make it possible to establish their connections with geophysical fields. Within the framework of the complex, it is convenient to combine heterogeneous data and integrate geophysical methods.
The article continues the series of works of the authors devoted to the deep structure of the White Sea region [21,22]. Its purpose is to model the velocity structure of the region's earth crust based on instrumental observation data using the Integro software complex. The objectives of the research are the construction of the velocity layers of the earth's crust, and the study of their connections with density heterogeneities and geophysical fields.
Density models were developed based on a reference velocity model of the crystalline crust, which has proven itself in a number of studies [23][24][25].
The primary data of seismic profiles on the structure of the earth's crust was refined within the framework of 2D density models and, using the 3D model, these results were extended to the earth's crust of the entire region. The relationships discovered during the modelling process are applicable in establishing the correspondence of structures at different depths, in the development of mineragenic criteria for mineral prospecting.

Initial Data
The geographical position of the White Sea region, and the structural elements of the sea and land area are presented in Figure 1.
Geophysical surveys in the White Sea Throat and Funnel by means of the CDP seismic reflection method, seismic refraction method, gravimetry, magnetometry, and seismoacoustic profiling were carried out by the OAO Marine Arctic Geological Exploration Expedition and FGUP Sevmorgeo, and in the adjacent land areas by Sevzapgeologii, Spetsgeofizika, "Geon Center", and "KSC RAS" [21]. In 2004, OAO MAGE performed the CDP seismic reflection method for a distance of 1700 km [26] on 11 profiles in the White Sea basin. The length of the active portion of the digital receiver was 3000 m, the number of seismic channels was 240, the distance between the centers of the groups of seismic channels was 12.5 m, the multiplicity of observations was 60, the length of the record was 6 s, and the sampling step was 2 ms. The source of elastic vibrations was a line of 10 Bolt-1900LL-XT-AT pneumatic sources with a total volume of 1500 cubic inches and an interval of excitation of elastic vibrations of 25 m. In addition, FGUP "Sevmorgeo" performed work by the deep seismic sounding (DSS) method on the 3-AR profile with bottom stations spaced 5-10 km apart. The length of the hodographs reached 300 km. The seismic sections used in the modelling process represented vertical sections of the lithosphere digitized at the appropriate scale with the profile to depths of 50-60 km, with boundaries separating areas with different velocities of P-wave propagation. The boundaries set the initial frame of the block structure, which was further refined taking into account the anomalous gravitational field. Constructed according to the data of gravimetric surveys [27-29], a digital map of the anomalous gravitational field in the Bouguer reduction at a scale of 1:1,000,000 ( Figure 2) was used to solve direct and inverse problems concerning the relationship of the gravitational field and the density distribution of rocks. The relief of the earth's surface was characterized by land elevation and sea depth marks.    Geophysical surveys in the White Sea Throat and Funnel by means of the CDP seismic reflection method, seismic refraction method, gravimetry, magnetometry, and seismoacoustic profiling were carried out by the OAO Marine Arctic Geological Exploration Expedition and FGUP Sevmorgeo, and in the adjacent land areas by Sevzapgeologii, Spetsgeofizika, "Geon Center", and "KSC RAS" [21]. In 2004, OAO MAGE performed the CDP seismic reflection method for a distance of 1700 km [26] on 11

Methods of Data Processing
Based on these marks, the values of the Bouguer anomaly, and the profile depths of the Moho boundary, a linear regression equation was derived, which was subsequently used in one of the variants for constructing the Moho surface (Section 3.3). Another regression equation expressed the relationship between the densities of rocks ρ (kg/m 3 ) in the earth's crust and the velocity V p (km/s) of longitudinal seismic waves: And formed the basis of the reference velocity model. The calculated rock densities reference five layers of the earth's crust model with the bedrock layers boundaries K1-K3, M, as shown in Table 1 [30,31]. The method of modelling using the Integro complex included the selection of an environmental model and its geometric framework, construction of 2D density models from seismic profiles, and the transition to a 3D density model of the region's earth crust [32]. The geometric framework and the environmental model were selected based on the block structure of the seismic profile and the values of the block densities calculated from the reference velocity model as initial approximations. The initial framework schematically repeats the profile structure in the form of deformed sections of density layers of the earth's crust broken up by series of the faults.
Direct and inverse problems are solved on two-dimensional and three-dimensional grids in Integro [33]. The direct problem was solved using 2D density models. The values of the selected structure blocks' density varied within the specified limits, achieving a minimum difference in the values of the calculated and the profile Bouger anomaly. If it was necessary to improve the optimization, the blocks were further divided, new heterogeneities were introduced, and the block boundaries were shifted. When solving inverse problems, the Integro complex uses the method of regularization, advanced algorithms for discrete Fourier transform, interpolation, and extrapolation of grid data [34][35][36]. The algorithms eliminate the edge effects caused by the lateral limits of the models and the lack of data periodicity, and are fast enough for computer processing of large arrays of geological and geophysical data. The 3D density model was built by solving the inverse problem of gravity prospecting [25]. The material of the mantle was considered homogeneous, and the density heterogeneities inherent in the earth's crust. The 3D model is convenient for constructing sections and calculating the contributions of individual elements to the anomalous gravitational field. In this work, this model is used to determine the spatial position and visualize the boundaries of the velocity layers of the earth's crust.
When constructing the boundary surfaces of velocity model layers, the coordinates of sites with corresponding longitudinal elastic wave values were marked on the seismic sections of the profiles. The obtained vector data were interpolated within 2D models and then transferred to a 3D model using B-spline multilevel interpolation [36]. The problem areas of the intersection of the layers' surfaces were corrected, after which the data interpolation was repeated. Figure 3 shows a 2D density model of a fragment of AGATE 3 geotraverse (Ust Pinega-White Sea), demonstrating the layered-block structure of the earth's crust [18]. The fragment is oriented across the main deep faults at the contact of the eastern part of the Kola and Karelian geoblocks, and extends across the Arkhangelsk Scarp, the Keret-Pinega Graben, the Poltin-Yelkib Horst, and the Leshukov Graben from south to north. The trend of the gravity field and the density distribution of the rocks in the 2D model reflect the boundaries of these rock units. The profile clearly shows the boundaries of the basement (depth 1-4 km, V p = 6.0-6.3 km/s) and M-discontinuity (depth 36-40 km, V p = 8.0-8.2 km/s). Manifested in geophysical fields, heterogeneities refer to the upper and middle crust. The lower crust with P-wave velocity of 6.8-7.1 km/s has a nearly uniform distribution of petrophysical parameters.  The boundaries of the layers of the reference velocity model in Figure 3 are highlighted in red. Such boundaries for all considered profiles define the frame of the 3D surface of the layer. Figure 4 shows a 3D density model of the White Sea region (a) and the position of the reference The boundaries of the layers of the reference velocity model in Figure 3 are highlighted in red. Such boundaries for all considered profiles define the frame of the 3D surface of the layer. Figure 4 shows a 3D density model of the White Sea region (a) and the position of the reference velocity model's layer boundaries in space (b). To simplify the picture, the positions of the layer boundaries are shown as in-depth contours obtained by the Isolines procedure of Integro in Figure 5.  As the diagram of the K1 surface depth (Figure 5a) shows, the sedimentary layer is present in some parts of the region and is mainly traced in the directions of the inflow of sedimentary material by the Severnaya Dvina and Mezen rivers. There is no sedimentary layer in Onega and at part of the Kandalaksha bay, near the Tersk Shore of the White Sea. Its maximum thickness of 6 km is observed at the confluence of the river Mezen in the White Sea and the White Sea Funnel.  As the diagram of the K1 surface depth (Figure 5a) shows, the sedimentary layer is present in some parts of the region and is mainly traced in the directions of the inflow of sedimentary material by the Severnaya Dvina and Mezen rivers. There is no sedimentary layer in Onega and at part of the Kandalaksha bay, near the Tersk Shore of the White Sea. Its maximum thickness of 6 km is observed at the confluence of the river Mezen in the White Sea and the White Sea Funnel. As the diagram of the K1 surface depth (Figure 5a) shows, the sedimentary layer is present in some parts of the region and is mainly traced in the directions of the inflow of sedimentary material by the Severnaya Dvina and Mezen rivers. There is no sedimentary layer in Onega and at part of the Kandalaksha bay, near the Tersk Shore of the White Sea. Its maximum thickness of 6 km is observed at the confluence of the river Mezen in the White Sea and the White Sea Funnel.

3D Density Model of the White Sea Region
The bedrock surface of the "granite-metamorphic layer K1-K2" (Figure 5b) is lowered in the southwestern and northeastern directions, in the Keret, Padun, and Leshukon grabens. In the central part, the depression is confined to the sedimentation basins of the Severnaya Dvina and Mezen rivers. The "granulite-basic layer K2-K3" (Figure 5c) varies in thickness from near-zero values to 30 km. A large area is occupied by the southwestern depression of its lower boundary. The Zimny Bereg Uplift is represented by a zone of differently directed gradients on its lower surface, and the depression is outlined in the White Sea Throat, which contrasts more in the topology of the M-discontinuity (Figure 5d).
On the M-surface, the depression in the White Sea Throat is surrounded by uplifts near the Mezen, the Tersk Shore of the White Sea, and Kanin Nos Peninsula. The maxima of M-depression are consistent with positive Bouger anomalies and the minima with negative anomalies (Figure 2).

Alternative Map of the M-Border
To check the simulation results, an attempt to obtain some of them independently was made. It is most convenient to carry out such a check on the example of M-discontinuity. Due to the density jump, this boundary is confidently fixed on seismic profiles, and its features are manifested in the anomalous gravity field in the Bouger reduction, overlying layers of the earth's crust and the relief of the earth's surface [37]. A linear regression equation: is calculated by connecting the depth Z (km) of the M-discontinuity on available seismic profiles, the value of the anomalous gravity field in the Bouguer reduction ∆G (mGl), and the elevation of the relief h (m) (Figure 6a,b). Using Equation (2), an alternative map of the depth of the M-discontinuity surface was constructed and the result is shown in Figure 6c. The map, with smoothed local features, in general correctly reflects the general plan of the structure of the M-discontinuity shown in Figure 5d.

Discussion
The well-defined crustal structure of a 3D model of the region in the form of alternating variably deep layers with differences ( Figure 4b) is consistent with generally accepted geological concepts.
The real structure of the earth's crust is complex. Instrumental measurements indicate the existence of high-speed interlayers in low-density layers, as, for example, in Figure 3c, or asthenospheric lenses at the boundary of the earth's crust and mantle. Nevertheless, the obtained 3D density structure of the earth's crust can be used to describe the deep structure of the region, clarify the structural boundaries, and develop criteria for identifying promising areas for the search for minerals.
The absence of the sedimentary layer in a part of the White Sea (Figure 5a) may be caused by the uplift of these parts associated with the Fennoscandian shield. A thick sedimentary layer in the White Sea Funnel hinders deep water exchange in the White and Barents Seas and is one of the causes of the low salinity of the White Sea, about 26 ppm, and associated ecological and biological effects.
The uplifts of the underlying surface of the "granite-metamorphic layer" fall on the structures of the Karelian and Kola megablocks, the Zimny Bereg Uplift, and the Onega peninsula ( Figure 5b). In this connection, the Zimny Bereg Uplift and the Onega Peninsula can be considered as the continuation of the structures of the Tersk coast of the White Sea and the Karelian megablock, respectively.
An example of the connection of the layer's form with geological processes is provided by the depression of the M-discontinuity, identified in the White Sea Throat. The results of modelling generalized empirically show that there is a subvertical crustal structure near the Zimny and Tersk Shores, the White Sea Throat, and the White Sea Funnel associated with kimberlite magmatism. Diatremes, sills, and explosion pipes are located in some parts of the Archean-Early Proterozoic crust with a mosaic type of geophysical anomaly and are confined to linear tectonic zones [12].
The depression of the M-discontinuity in the White Sea Throat, reflecting earth crust-mantle interaction, is a criterion for diamond potential [38]. The similarity of its local extrema (Figure 5d) indicates the interpenetration, mixing, and transformation of deep-seated matter in the conjugation zone of the uplifting Fennoscandian Shield and the Russian Plate. Geodynamically, this process is contrary to tectonic plate subduction. The Zimny Bereg Diamondiferous Province is associated with an old crystalline basement scarp, which is constantly uplifting. The Zimny Bereg kimberlite field is located in a NW-trending Riphean aulacogen, in the Kola Craton-Mezen syneclise conjugation zone. The gravity field of the aulacogen passes from weakly negative to positive values. The magnetic field values are elevated. Gravity and magnetic anomalies in lower crustal layers, affected by low heat flow, indicate a deep source of magma generation. Regional faults and areas of intersection of tectonic zones display deformations, fastenings, the steep bends of the contour lines of gravity and magnetic fields, and the geochemical aureoles of Ni, Mn, Cr, Co, and Zn. Ultramafic explosion pipes, occurring in the crystalline basement with high electrical resistance, exhibit elevated conductivity [12,19].

Conclusions
Two-and three-dimensional density models of the White Sea region's crust were developed based on the results of well-known seismic and gravimetric instrumental studies using the Integro software complex. In the 3D density model, the positions of the boundaries of the speed layers of the earth's crust of the region were determined.
Boundary surface relief schemes allow the local features of layer's boundaries to be geographically referenced, features to be associated with gravitational field anomalies, and the power of Earth's crust layers to be estimated.
The proposed interpretation of the simulation results explains the absence of a sedimentary cover on part of the water area and its large capacities in the deep-sea areas by raising of the Fennoscandian shield. Connections within the framework of the granite-metamorphic layer of the Winter Bank Uplift and the Onega Peninsula with the structures of the Terek Coast and the Karelian Megablock, respectively, are indicated, in addition to large capacities of the granulite-basite layer in the area of Winter Bank lifting by the intensity of its formation, and the depression of the M boundary by the existence of a subvertical structure associated with manifestations of kimberlite magmatism.