ARCTIC REGION AS THE BASE FOR GEODYNAMIC INTERPRETATION

Cluster analysis is applied for computing stable combinations of geological and geophysical parameters, and areas with such combinations are interpreted as regions that differ in structural and geodynamic features. The shelf areas are distinguished by specific sets and patterns of parameters, including sedimentary cover thickness, tectonic heterogeneity of the basement, heat flow, anomalous magnetic field, and gravity anomalies that reflect the topography of the crust–upper mantle boundary. In the deep oceanic areas, S-wave velocity variations show abnormally ‘cold’ blocks, while the average heat flow values are high. This combination of parameters is typical of transform zones at the junction of the Atlantic and Arctic segments. Superimposed thermal domes are located symmetrically with respect to the axis of the mid-oceanic ridges (MOR). Such domes may occur on the continents located close to MOR. Similar indicators can be revealed along the transition zone to the north of the East Siberian Sea.


MULTIDIMENSIONAL DATA ANALYSIS, AND AN APPROACH TO SOLVING THE PROBLEM
Today the geological and geophysical information on the Arctic floor structure ( Fig. 1) is abundant, and attempts to compare and analyze the available multidimensional data files by conventional visual and correlation methods do not seem practical. Despite the fact that the Arctic shelf area has been unevenly covered by the studies, the authors attempt at generalization of the source data and apply the cluster analysis to calculate and classify combinations of geodynamic parameters.
Cluster analysis results can provide the basis for mapping the studied territory by geodynamic settings and interpreting the physical meanings of the identified types. So, the objectives are to definite the terms of 'geodynamics' and 'geodynamic setting', select geological and geophysical parameters to describe a 'geodynamic setting' (providing a 'spatially uniform' description whenever possible), select a computation technique for identification of geodynamic setting types, implement a computational algorithm and describe its specific features, construct a map showing patterns of geodynamic setting types, and interpret the physical and geodynamic meaning of the results.

THE TASK OF GEODYNAMIC ZONING, AND EARLY WORKS
By comparing the concepts proposed by researchers holding diametrically opposite viewpoints [Belousov, 1975;Pavlenkova, 1987;Zonenshain, Kuz'min, 1993], one can find that there is the single 'core' based on the definition of 'dynamics' as a discipline, which is accepted in physics. "Mechanics studies the simplest form of motion of matter -mechanical motion, i.e. changes of the relative positions of bodies or their parts in space in the course of time. Bodies are macroscopic systems consisting of a very large number of molecules and atoms, and the sizes of these systems are many times larger than the intermolecular distances. Kinematics studies mechanical motion of bodies without reference to the interaction between the bodies. Dynamics considers the effect of the interaction between the bodies on their mechanical motion" [Yavorsky, Detlaf, 1974, p. 13-14]. 'Interaction' envisages analysis of forces and energy sources. Thus, in our study, 'geodynamics' is considered as the science that studies the interaction between geological objects in the course of time.
An essential issue in problem solving in geodynamics is developing approaches to define parameters that can describe complex properties of geological bodies for further quantitative analyses. The definition of 'geodynamics' in [Khain, Lomize, 1995, p. 5] is very close to the above citation. Objects in terms of geodynamics are beyond measure more complicated than objects in classical physics, which makes geodynamics a unique discipline. Due to the high complexity of the objects, it is challenging to develop adequate and effective quantitative models that can properly describe the geodynamical processes.
With regard to the above definitions, it is reasonable to distinguish three main groups of parameters describing a geodynamic object: (1) parameters describing geometry and physical properties; (2) parameters describing forces and energy release; (3) parameters characterizing motions resulting from impacts of forces on the object and the energy release in the object.
So, the geodynamic zonation objective is to search for different stable combinations of parameters describing geodynamic objects and analyze the patterns of such combinations in space.
The Arctic region geodynamic is described in many publications, including papers based on modern seismic tomography data of higher resolution than in the late 1990s and early 2000s [Jakovlev et al., 2012;Koulakov et al., 2013]. Studies of the East European Platform and the Caucasus reported in [Reisner, Reisner, 1987, 1990 pioneered in applying the formal quantitative approach to solution of geodynamical problems using several parameters. The cluster analysis results for Western Eurasia were published in [Ioganson, Boltyshev, 2000]. The first publication on the selection of parameters and calculations of cluster combinations for the deep areas of the Atlantic Ocean was [Sokolov et al., 2008] (although measurements of vertical movement velocity, one of the most important parameters, were lacking).

THEIR COORDINATES
Selected parameters should be defined, where possible, to the same detail across the study region to ensure adequate comparison and evaluation of different areas identified in the study region. However, in practice, it is not feasible due the non-uniform coverage of the areas by studies and the incompleteness of the available data sets. For computing, a full array of parameters should be assigned to each cell in the grid (with the uniform spacing along X and Y axes), and if measurements are sparse, values interpolated over long distances should be used. Heat flow is the parameter that fits best for geodynamic calculations, although it is least studied in the Arctic. The density of seismic tomography data is strongly dependent on the density of events and the number of stations, but the stations are rare and widely scattered in the Arctic. The selection is also limited by the need to use P-and S-wave data of similar detail. A well-detailed database, including data from 3D surveys, is available for the sedimentary cover in the shelf zones, while the database for the deep Arctic basin is incomparably smaller. In this regard, we have to refer to the data set which level of detail for shelf areas is reduced to ensure comparability to deep oceanic areas. The most detailed data set is available for the gravity field from satellite altimetry, although the signal-noise ratio decreases abruptly for wavelengths less than 25 km [Sandwell, Smith, 2009]. This predetermines an approximate minimum acceptable step of the grid. Using larger cells leads to gradual delimitation of the shelf-deep ocean transition zone, but this zone is important, and it is thus reasonable to put a limit of 25 km, while realizing the redundancy of this quantization step for many of the less detailed parameters. Anyway, this is acceptable. Besides, this step makes up a grid that can be potentially supplemented as and when more detailed information layers become available. The 25-km grid provides for constructing a 1:10000000 map.
The selected parameters should describe the three groups of properties mentioned in Section 2.1. It is easy to select parameters describing structural features of the lithosphere (Group 1) (see Section 3). Parameters defining energy release (Group 2) can also be easily selected, although there is a problem of irregular heat flow measurements.
The major issue is how to select parameters describing the resultant motion (Group 3). For the land areas, velocities of vertical movements recorded by repeated geodetic measurements and GPS data can be referred to. Such measurements are not available for the sea floor, and a regular observation and measurement grid can hardly be obtained in the near future. Therefore, to include parameters of Group 3 into calculations, there is a need to use so-called 'surrogate'' parameters that can indirectly reflect unmeasurable values or partial quantities measured in the ocean, or those reflecting a combination of many effects including those to be processed. This approach is the only way to reflect the necessary information in the absence of detailed data. In this paper, only the recent state of the lithosphere is discussed. Values of the parameters in the three groups are obtained by instrumental measurements. Particularly, it should be noted that among the analyzed parameters, there is no 'pure' parameter that would belong clearly and only to one of the three groups.
The recent progress in the Arctic shelf zone studies is significant, but our main objective is to study the deep oceanic areas of the Arctic and find a possibility for comparing them with the Atlantic, and the analysis depends on the current coverage of the areas by the studies. It is thus reasonable to limit the level of details in the shelf database with regard to the selected grid and interpolate it to the rare deep ocean measurements that are currently available. In the future, results will have to be revisited and updated with account of newly obtained measurements. As the Arctic shelf territory is wide, the cluster analysis of this region (as opposed to the Atlantic) [Sokolov et al., 2008] covers a wide range of morphological structures associated with the continent-ocean transition zone, with account of the fact that the shelf is a continental structure. Thus, the data on deep ocean basins are analyzed, and the 'land-sea' data are linked in our study. The most promising development of such studies can be cluster analysis of data on West Siberia and the territories towards the Kara Sea, taking into account the complete database availability and the importance of these territories.

APPROACH TO SELECTION OF DATA PROCESSING METHODS
Among the multivariate statistical classification techniques, there are three techniques that have been successfully applied to classify the geological and geophysical -discriminant, factor and cluster analysis, that differ in specific features. Discriminant analysis refers to the known a priori stable type and classifies objects by comparing the values of parameters in the specified set against the 'standard'. Factor analysis assumes that the available data set reflects a combination of the effect of two or more processes, each contributing to values of all the parameters. Cluster analysis identifies stable combinations of parameters which are not detectable by visual analysis of maps. This technique is the most adequate at this stage of our studies, judging from results of its application in the studies of other regions (see para. 2.1).

DATA OVERVIEW
The Arctic Ocean shelf and floor are objects of geodynamic zoning (Fig. 2). Our calculations are based on the following parameters: oceanic floor relief, sedimentary cover thickness, surface Love wave tomography, Bouguer anomaly, isostasy, heat flow, S-and P-wave seismic tomography, total seismic moment, and magnetic field anomalies. To calculate cluster combinations, each cell of the grid is assigned an array of ten parameters defining a point in the multidimensional space (after normalization and centering control). The difficulty is that a space projection of every axis to each other is zero, i.e. the parameters plotted along the axes are linearly independent. For the available geological and geophysical parameters, this rule is not satisfied.
Correlations between the parameters are not zero, yet not close to 1 or -1, so all the parameters are, to some extent, linearly dependent. Such an issue is due to the fact that one and the same substance of the lithosphere is the source of different fields (a direct problem), and a cross-correlation between these fields is most probable. In our study, there is an inverse problem (patterns of sources of the known fields are unknown), so our calculations are based on what is available -sets of linearly dependent field measurements.
The selected parameters are briefly described below. Parameter values in each cell (25×25 km; polar stereographic projection) are used as components of multidimensional arrays for statistical processing.

OCEANIC FLOOR RELIEF
Oceanic floor relief is the first and one of the most important parameters describing the top layer of the crust and lithosphere (see Fig. 2). In our study, we use data from [IBCAO, 2008], smooth the values by low-frequency filtering and recalculate for the cell (25×25 km). The oceanic floor relief reflects impacts of many processes, including magmatism, deformation, sedimentation etc. In our classification, the relief is a directly measurable parameter in Group 1 (geometry). Its quantity is an indirect reflection of crustal block movements under the impact of forces applied to the crust (Group 3 parameter). Precise motion monitoring data (similar to GPS measurements on the continents) are not available for the oceanic floor. Anyway, the relief data processing ensures that results of the movements are indirectly (although still inadequately) taken into account.

SEDIMENTARY COVER THICKNESS
The sedimentary cover thickness is reconstructed for the Arctic region and the adjacent land ( Fig. 3) from the data published in [Laske, Masters, 1997]. In our study, the sedimentary cover thickness values are consolidated and averaged for a grid of 30 arc minutes in order to introduce relevant corrections in tomographic models. The values are adjusted to the working detail of this work. This isopach map is the only one showing the entire regions in a more or less uniform detail.
Taking into account the sedimentary cover thickness is important for several reasons. The main reason is that the ocean periphery is the zone of intense sedimentation (due to sediments wash-out from the continent), and the isostatic equilibrium between the crustal blocks and the viscous mantle is disturbed because of the increased load on the latter. This leads to processes aimed at restoration of the balance by relevant vertical movements so that the medium could achieve the equilibrium state and the disturbance would be smoothed  [IBCAO, 2008], and the region contours (black line) for multivariate statistical analysis.
out. Another reason for including the sedimentary cover thickness in the calculations is that densities of the bottom and crystalline base are significantly contrasting, and this surface must be taken into account when describing the properties by parameters in Group 1 (geometry of the object; in this case, a combination of the crust and the upper mantle).

SURFACE LOVE SURFACE WAVE TOMOGRAPHY
The surface Love wave tomography (Fig. 4) is based on the data from [Larson et al., 1999]. This parameter is 'surrogate' to describe the object's geometry as it cannot directly reflect the behavior of an effective bottom of the crust-mantle layer (direct depth measurements are not available for the mantle margins across the entire Arctic), but still provides some indirect information. Phase velocities of surface waves (i.e. waves propagating in the effective surface layer) depend on the layer thickness. The larger is the thickness, the slower are velocities of wave propagation, and vice versa, the smaller is the thickness, the higher are the velocities. Thus, a percentage deviation of the phase velocity from the average value reflects a relative variation of the surface layer, which is proportional to the required parameter (the depth of the crust bottom or the intra-mantle boundary). In the model chosen for our calculations. the wave period is 35 seconds, and displacements along the wave front penetrate to a shallow depth that roughly corresponds to the top of the lithospheric layer. Deeper layers are subject to wave motions with large periods. In Figure 4, clearly identified are continental areas and zones with lower velocities in the vicinity of ancient shields, which are distinguished from the oceanic zones with the thin and highvelocity lithosphere. A significant correlation is noted between the continental areas and the negative tomography field. A principal coincidence is revealed by comparing the anomaly's zero value position and the shelf margin shape. Moreover, this suggest two interesting exceptions to the general rule. First, the 'continental' values of the anomaly are revealed considerably far (about 450 km) into the oceanic area between the Lomonosov and Mendeleev ridges. Since the horizontal calculation accuracy is 200 km, it is impossible to specify exactly from the shape of this anomaly to what extent the continental crust block may be traced. Nonetheless, the presence of the continental fragments is evident in this area. The second exception is the presence of the 'oceanic' value of the anomaly underneath the Svalbard archipelago and the western Franz Josef Land (Fig. 4). It is known that this region is subject to intensive processes, including positive vertical and horizontal movements, Quaternary volcanism, formation of submeridional magnetic anomalies, increased heat flow. This may suggest a new stage in the  [Laske, Masters, 1997].
tectonic development of the Eurasian periphery in this area. In view of the above, the selected parameter can fit as a 'surrogate' to describe the geometry of the crust bottom and the shallow mantle margin.
Group and phase velocities of surface Love waves and S-wave velocities are not independent parameters -the velocities depend on the properties of the matter in the surface layer (0-100 km), and such properties are given by elastic module values for the same volume. Nonetheless, anomalous fields of these parameters are revealed in various patterns in different deep oceanic areas. As noted above, purely linearly independent parameters are absent in the calculations. So, seismotomographic data of all types are used, despite the correlation. Variable surface wave velocities in different periods cannot be precisely linked to variations in the surface shapes within the lithospheric layer. However, mapping of short-period phase velocities (in particular, 35 s phase) shows an evident relationship between surface wave anomalies and crustal thickness variations. Some exceptions are mainly associated with active magmatism (Iceland) and non-standard tectonic settings (Quaternary volcanism in the Svalbard archipelago shelf area). Anyway, in general, this parameters is geometrically proportional to boundaries inside the lithosphere.

BOUGUER ANOMALIES
Bouguer anomalies are calculated from AGP data, Arctic Gravity Project [Forsberg, Kenyon, 2005] and the topography data from [IBCAO, 2008] with account of average densities of the oceanic crust and the land (2.8 g/cm 3 , and 2.67 g/cm 3 , respectively) and 166-km radius integration (Fig. 5). AGP gravitational anomalies are values in free air. For the ocean areas, this means that about 80 % of the anomalous field variability is proportional to the most contrasting density boundary, i.e. the oceanic floor relief known from echo sounding surveys. In the Bouguer anomaly calculations, an impact of the floor relief on the anomalous field is eliminated by 'adding' a mass to the aqueous layer to reach an average value estimated for the crust. After such a procedure, the residual field variability reflects mainly the depth of the density-contrasting crust-mantle boundary and lateral density inhomogeneities in the crust and mantle. Such variations may be insignificant in the deep oceanic basins and quite significant in areas of serpentinization of the upper mantle rocks. In the absence of seismic survey data for the deep oceanic areas, impact of the variations cannot be reliably distinguished from variations in the crust bottom depth. At the same time, lateral inhomogeneities may have  [Larson et al., 1999]. The shelf edge is shown by the dashed line.
large spatial sizes and reflect heated zones in the lithospheric areas containing magma chambers and considerable partial melting. Such zones are marked by active magmatism and correspondingly increased thickness of the crust, which loads the viscous mantle and leads to deeper positions of the Moho. In view of the above, the Bouguer anomalies are generally proportional to the crust-mantle boundary depth, and the lower is the anomaly value, the deeper is the boundary.
In considerations of the exceptions requiring a thermal correction of the anomalous field, there is a need to take into account either heat flow or another parameter (e.g. S-wave tomography) in order to reflect the state of heating. However, these parameters are independently involved in our analysis (see below), so for the geodynamic analysis in this study, it is impractical to calculate the Bouguer anomalies with account of a thermal correction.
The Bouguer anomaly is are viewed as a parameter in Group 1 describing the geometry of the crust-mantle layer, i.e. the interior boundary separating 'dense layers' of the crust and upper mantle, as well as mass variations along the surface layer. These properties in the Bouguer anomaly are present as combined contributions that cannot be reliably separated. Thus, the Bouguer anomaly is a 'surrogate' parameter describing the geometry and a direct parameter describing the mass variations. Its geometric component is more sig-nificant as the density variations in the mantle and the lower crust layers cause less variations in the field than positions of the boundary. The crust thickness can be directly calculated as a linear function of the Bouguer anomaly (according to [Deminitskaya, 1967, p. 27] which results from approximations by comparison points to real DSS data. It is a fully linearly dependent parameter having the Bouguer anomalous field morphology. In our calculations, it is not used because the initial field is preferable when the calibrated DSS data coverage is sparse.

ISOSTATIC ANOMALIES
Isostatic anomalies are calculated from the data on the Bouguer anomalies and the topography data from [IBCAO, 2008] with account of average densities of the oceanic crust and the land (2.8 g/cm 3 , and 2.67 g/cm 3 , respectively) and 166-km radius integration in the Airy model with the reference surface depth of 33 km (Fig. 6). Long-wavelength components of more than 900 km are removed from the anomalous field as they reflect sublithospheric inhomogeneities, and their effect in isostasy conceals the processes taking place in the Earth's upper shell. After eliminating the anomalous field variability associated with the upper crustal boundary and calculating the Bouguer anomalies, isostatic anomalies are calculated to eliminate a hypo-  [Forsberg, Kenyon, 2005] and the topography data from [IBCAO, 2008], average density values of the oceanic crust and the land (2.8 and 2.67 g/cm 3 , respectively), and a 166-km integration radius.
thetical field variability associated with variations in the compensation surface topography due to variable thicknesses of the crustal blocks 'floating' on the viscous mantle surface. We assume that in case of isostatic equilibrium, the compensation surface position is a topography factor: where H is depth of the compensation surface, T is reference depth, σк is crust density, σв is water density, and σм is mantle density. Relevant correction factors are estimated accordingly to eliminate the effect of a hypothetical surface topography (in the same way as the floor relief impact is eliminated).
The residual field reflects isostatic anomalies which positive values indicate the presence of an excessive mass over the compensation surface, while the negative values show the lack of a mass. The excessive mass leads the crustal block slinking at the given point; if the mass is lacking, the block is uplifted together with the portion of the mantle. This is true only for cases envisaging a complete action that disturbed the isostatic equilibrium of the crustal blocks system. If the action (such as, for example, thrusting of one block onto another) is not completely fulfilled, both the excessive mass (i.e. positive isostatic anomalies) and positive vertical movements of the crust are present. So, in isostatic anomaly interpretations, there is an ambiguity that can be eliminated if more information is available and the general tectonics of the region is known. In terms of geodynamics, this parameter can directly describe variations in density properties of the crust, intensity of energy release in the crust, and stresses (isostasy gradient module) due to transition of the medium from the disturbed state to equilibrium. This parameter is also a 'surrogate' one for describing vertical movement of the crustal blocks which result from the energy release. These properties are reflected as an inseparable combination in the isostatic anomaly field.
In isostatic anomaly calculations, the gravity effect from the mantle surface is estimated from the floor relief data in proportion to the ratio of densities (Airi model). It is then subtracted from the Bouguer anomalies calculated from the same topography data. In this case, however, a full correlation between them is absent due to the fact that in 3D calculations, the topography effect in the Bouguer anomalies (for instance, in  [IBCAO, 2008] for average density values of the oceanic crust, the land and the mantle (2.8, 2.67 and 3.3 g/cm 3 , respectively), and a 166-km integration radius in the Airy model with the reference surface depth of 33 km.
Рис. 6. Изостатические аномалии, рассчитанные авторами по данным аномалии Буге и рельефу [IBCAO, 2008] для средней плотности коры океана 2.8 г/см 3 , плотности суши 2.67 г/см 3 и плотности мантии 3.3 г/см 3 при интегрировании с радиусом 166 км по модели Эйри и глубине поверхности приведения 33 км. the aqueous layer approximation by prisms in the grid nodes) is calculated as the integral with weights equal to a distance from the prism bases. By calculating the integral effect of the mantle surface reconstructed for the same topography, but approximated by prisms with depths relative to the reference surface (33 km), a smoother adjustment field is obtained as the surface is more distant from the sea level, and at the observation point there is no contrast weight contract that is present in the aqueous layer close to the observation point. Besides, the analysis of the Bouguer anomalies and isostasy at sublatitudinal profiles in the Atlantic [Sokolov, 2015] shows that these fields are positively correlated in the abyssal basins with the practically undisturbed isostatic equilibrium, and the correlation is lacking in the rift zone and at the ridge shoulders up to 300 km wherein the isostatic equilibrium is disturbed by rifting.

HEAT FLOW
The heat flow reconstruction (Fig. 7) is based on the data published in [Hasterok, 2011]. Heat flow data are very randomly available for the Arctic aquatic area, but this parameter is absolutely essential for our calculations, so we have to use whatever numerical values are published. Since clouds of values assigned to each cell are uneven and highly scattered, the grid is estimated by kriging, followed by smoothing high-frequency components to the level of other parameters in order to minimize the effect of the uneven density of measurements. The resultant map in Figure 7 reflects the uneven scatter of the data. Heat flow is a direct parameter in Group 2 (energy release). In view of the uneven coverage of the Arctic by studies, it is obvious that areas with the dense observation network can be reliably classified by this parameter, while almost no impact of this parameter is revealed in poorly studied areas.
To study the crust and upper mantle geodynamics, we select the topmost section of NGRAND model (0-100 km) that was calculated by its authors for 2×2° blocks and represented by spherical harmonics up to order 31. The tomographic matrix is recalculated for a 25×25 km grid. The values represent percentage deviations of the transverse wave velocities from the average value estimated for the layer. Estimating this parameter is sensitive to heated areas and zones with considerable partial melting. The parameter perfectly shows the presence of pluming (usually accompanied by magmatism) and ridge zones. In these zones, tomography values are negative (for deep plumes, -3.5 % and below)  [Hasterok, 2011]. The mid-oceanic ridge axis is shown by the black line.
as seismic wave velocities are decreased in the heated and less viscous medium. So, this parameters is a combined reflection of energy release effects (in the heated medium) and the medium's geometry (zones with active magmatism and increased thickness of the crust), and this combination can hardly be split. In the two groups of parameters, it is a 'surrogate' reflecting the group properties only indirectly.

P-WAVE TOMOGRAPHY
P-wave tomography ( Fig. 9) is based on the data published in Becker, Boschi, 2002]. To study the crust and upper mantle geodynamics, we select the topmost section of HWE97p model (0-100 km) that was calculated by its authors for 2×2° blocks and represented by spherical harmonics up to order 31. The tomographic matrix is recalculated for a 25×25 km grid. Similar to the case of transverse waves, longitudinal waves should reflect the thermal state of the interior, and P-wave tomography should correlate with S-wave tomography. However, in practice, this is not the case. According to [Becker, Boschi, 2002], the Sand P-wave models are better correlated towards the middle part of the mantle (over 1000 km), which demonstrates that the parameter variability causes are similar. In the upper mantle and at the surface, the Sand P-wave models are considerably different. The behavior of the S-models can be explained in a consistent way, but the pattern of values in the P-models can be explained only if other options are suggested for interpreting the velocity variation sources. In our opinion, the variations may be caused by stresses and/or fracturing in the lithosphere, which result in a quire specific pattern of maximum and minimum values of this parameter, as shown in the map (Fig. 9). The minimum values are concentrated along collision zones of the Earth. The maximum values are revealed rearwards of the collision zones, and microfractures in such rear areas are not coincident in plan with directions of forces that form the collision zones. It is probable that collision is the factor of 'hindrance' for longitudinal waves. This parameter reflects the state of stresses of the medium and the corresponding fault system that releases the stresses. It can thus be viewed as a combination of parameters (energy release, and forces impact results) and a 'surrogate' in both Groups 2 and 3. A consistent regional geodynamic interpretation of this parameter is lacking, so considering it in the context of other parameters is even more interesting.
P-wave velocity variations in the lithospheric layer of the Arctic are shown in the modern models [Jakovlev et al., 2012] in much more detail than in HWE97p model. It is desirable to use such models for multivariate statistic calculations in combination with similar detailed models of other tomographic parameters, such as variations of S-wave and surface Love wave velocities, but this is unrealistic because of the lack thereof. It  [Grand et al., 1997;Becker, Boschi, 2002]. The shelf edge is shown by the dashed line.
is noteworthy that the datasets similar to those used in our study have been already used in the calculations for the Atlantic lithosphere [Sokolov et al., 2008], and the calculation area of the northern Mid-Atlantic ridge is overlapping with that of the Arctic region. Therefore, it becomes possible to conduct the cluster classification for subregions of the uniform Atlantic-Arctic global ocean segment.

TOTAL SEISMIC MOMENT DENSITY
This parameter (Fig. 10) is calculated as the total energy released by earthquakes from the dataset on seismic events (M>3) in the 0-100 km thick layer northward of 60° [ANSS, 09.12.2010]. The calculation method is described in [Boldyrev, 1998]. Energy amounts released by the events in the given cell are summed up as follows: М=(10 (17.1+1.3·(Mag-5)) )/10 +13 [joule·10 +13 ], where M is total moment, Mag is magnitude (Richter scale). After calculating the total moment, the seismic moment density per square kilometer is estimated for each cell with regard to its square area. The final dimension of the value shown in the map (Fig. 10) is [j/ km 2 ]·10 +13 .
Values of this parameter are very unevenly scattered in the study region. Moreover, the energy release along MOR exceeds 5 % of the global seismic energy. So, variability of this parameter is mainly manifested outside of the study region. In this respect, it differs from other parameters which values close to absolute minimum and maximum are revealed within the study region. It is noteworthy that any limits to registration of seismic events by distance are not specified for the selected magnitude threshold. It is thus acceptable that the seismic moment density values are uniformly scattered across the entire region, although equal to zero in the major part of the region. This parameter is a direct one in Group 2 (energy release).

ANOMALOUS MAGNETIC FIELD
This parameter is derived by EMAG2 data processing [Maus et al., 2009] (Fig. 11). Maps showing anomalies of the complete array, vertical component and full gradient module of the anomalous field are available in [Maus et al., 2009]. In our calculations, the latter parameter (analytical signal) is used due to its advantage -the lack of any alternating-sign field resulting from the alternating directions of the magnetizing field. It provides the most proper reflection of the crust properties, which facilitates the data interpreting. This parameter is proportional to the concentration of magnetically active minerals in the lithosphere and reflects the factors controlling the concentration variability (Curie isotherm depth, serpentinization zones, intense Fig. 9. P-wave tomography the 0-100 km layer, based on the data from Becker, Boschi, 2002]. The shelf edge is shown by the dashed line. Рис. 9. Томография по Р-волнам для слоя от 0 до 100 км по данным Becker, Boschi, 2002]. Пунктирная линия -бровка шельфа. magmatism zones that differ in composition from the surrounding areas, etc.). The anomalous magnetic field (AMP) is a 'surrogate' parameter in Group 2 (energy release) and, to some extent, in Group 1 (geometry of the deep boundaries).

CLUSTER ANALYSIS (GENERAL CONCEPTS)
Cluster analysis is a multivariate statistical classification technique that places measurements into groups (i.e. stable combinations of parameters in a multidimensional space), defines the geometry of such groups, estimates distances between group centers, and defines spatial limits for each group. An initial set of points in a multidimensional space (according to the number of parameters used for classification, i.e. 10-dimensional space in our study) is divided into clusters. Points in a given cluster tend to be similar to each other in some sense, and points in different clusters tend to be dissimilar. In our analysis, an object is a 25×25 km cell, and values of ten parameters are assigned to each cell. Generally, a cluster is viewed as a group of objects (in this case, the crust and upper mantle) which has the property of density, i.e. a compact concentration of the parameters in the specified area of space; the density of the objects (i.e. the simila-rity of properties) is higher within the cluster than outside it; each cluster has its center and dispersion (i.e. effective array); its shape is a hypersphere; a cluster is separable from other clusters. In fact, this definition of the technique is not precise, but clearly defines its capacities and objectives, most of which are intuitively comprehensible.
In our study, calculations are performed by STATISTICA software after loading the prepared data. This means that the authors are not aware of the details of the algorithms implemented in this software package and know only the general classification method controlled by parameters listed in the user menu. The main parameter is the number of clusters, N. Our intention is to break down the entire set of objects into this number of clusters. The procedure for selecting the optimal number of clusters is described below in para. 4.2.
The source data are standardized parameters (see Section 3) for each cell (see para. 6.1). In the table, the columns show values of each of the 10 parameters, and the rows correspond to the cells. In the multidimensional space, a matrix of distances between each pair of objects is calculated, the number of clusters (N) is set, and the algorithm begins to break down the entire set into clusters. In general, the procedure is as follows: the radius is set to exceed the size of the entire cloud of the objects so that the radius can reach from one object to another; the algorithm begins to decrease the radius until the cloud is broken down into separate dense group; and the access to each other with the current radius is no longer possible.
The above-described computerized technique is applied for zoning of the multidimensional space characterized by a high concentration of points. It is very much similar to visual and manual methods used for zoning the territories in conventional 3D and 2D studies. In terms of physics, its meaning is very simplewith its own computer code, it is not complicated by any additional data processing parameters, such as weights of factors. The technique is called k-means clustering, and it is unknown which specific calculation algorithm is implemented in STATISTICA software. It should be noted that the relevant software module has been successfully applied in our pilot study, so there is no need to use any other software packages or develop any proprietary modules.

APPROACH TO DEFINING THE RESULT-ACHIEVEMENT CRITERIA
As briefly described above, the objective is to classify all the objects into stable and clearly separable statistical groups, which number, N should be as large as possible, and each of such groups should be characterized by its own combination of all the selected parameters. It is obvious that groups with extreme values of any parameter are the first to be detected. The division by less pronounced variations in parameter values is implemented only after the groups with maximum values (and values in the main range of variability) of each parameter are identified. At this stage, it is important to fix the moment when the principal division of the specified areas by statistically distinguishable average values is stopped, and the 'forced' division is commenced to identify clusters that differ from each other by a small value commensurable to dispersion or an instrumental error of a parameter within a separate area. At this moment, linear estimations of the medium heterogeneity are complete, and the procedure goes on to analyze the scattered heterogeneity. In such a case, the geodynamic interpretation of individual clusters seems to become pointless, so the analysis should be stopped at the current value of N, and a variety related to the scattered heterogeneity should be assessed in statistical characteristics such as high-order moments, ensuring the uniformity across the entire area. Another result-achievement criterion is physical and geological substantiation of various parameters in each cluster. So, the problem of geodynamic zoning is considered solved when the specified characteristics sets are evaluated for each parameter in each cluster.

GEODYNAMIC CLASSIFICATION ALGORITHM
The procedure of data preparation for clustered classification is as follows: all the parameters (see . The anomalous magnetic field (AMF) gradient module according to EMAG2 data from [Maus et al., 2009]. The shelf edge is shown by the dashed line.
tion 3) are arrayed in matrixes which spatial dimensions are identical; values of the parameters are standardized as required for using the distance calculation algorithm (dimension of parameters should be similar); the data are tabulated; and the table is loaded for computer processing.
The next procedure is testing the classification with small N values -the algorithm is applied to conduct the step-by-step primary classification of the analyzed space into clusters with the obvious geological interpretation. Starting from N=2, the oceanic and continental (shelf) areas are differentiated. In the next step (N=3), the shelf area is sectored by values of the sedimentary cover thickness. Transition zones are consistently identified in the next steps up to N=5. In steps N>5, in addition to trivial solutions, it becomes possible to detect the settings that are not clearly detectable by the visual analysis. For instance, the shelf areas and basins are differentiated by heat flow values and tomographic parameters. In steps from N=8 to N=10, zones are differentiated by the magnetic field and Bouguer anomalies. The shelf areas are detected prior to the deep oceanic areas in the Arctic, which indicates that the shelf areas are less heterogenic in view of the selected parameters. Steps from N=11 to N=14 ensure the final stable differentiation of the deep ocean zones.
In steps from N>20, the spatially largest cluster of deep oceanic basins is 'broken down' into small clusters that are randomly scattered in space and differ from each other by amounts comparable to the average dispersion of the parameters in the standardized medium. At some point, the number of clusters is sharply increasing across all the Arctic areas, which profiles are concentrated in the field of zero scattering and do not deviate from zero by any significant values (Fig. 12). This shows that the physically justified limit of the cluster classification is reached for the available data set, and any further increase in N with its asymptotic approximation to the number of objects (or to infinity, depending on the detail of partitioning in the analyzed territory) will not improve the solution of the classification problem.

CLUSTERS OF GEOPHYSICAL PARAMETERS
The Arctic statistical clustering model is based on the parameters described in Section 3, estimations obtained by the described technique (Section 4), and the solution selection algorithm (Section 5). It comprises 14 stable combinations of the parameters listed in Table 1. The key point of the applied technique is evaluation of distances in the multidimensional space, and the calculations are done for standardized parameters given in similar measurement units with zero mean dispersion and zero unit dispersion. Among two options for calculating the normalizing statistical moments -only for the study region or for the entire Earth, the first option is selected because all the data sets used in this study contain extreme values close to the absolute minimum and maximum values recorded for various regions of the Earth, except for the total seismic moment which average values in the obtained clusters are by 4-5 orders lower than the maximum values estimated for island-arc zones of the Pacific Ocean. Nonetheless, this parameter is also normalized to the regional value.
In Figure 12, estimated central parameter values for the clusters are plotted in dimensionless coordinates. The profiles show that the cluster combinations contain values of each parameter, which reflect almost the entire main range of values, ±σ (equal to 1). All the main values are involved in one or another stable combination, i.e. cluster. The exception is the seismic moment due to its extremely uneven scatter and the presence of a sharp extremum, as well as asymmetrical tomography by various waves because the study region is predominantly oceanic.
At the final stage of the computerized data processing, the parameters combination patterns are mapped. Each cell (25×25 km) is assigned its own cluster number according to the calculation results, and the cell is painted in a color assigned to this number. The result of this procedure is the Arctic geodynamic zoning map (Fig. 13). This map and Table 1 are the main result of our study. Further discussions refer mainly to the map.
Estimated square areas of all the clusters are given in Table 2. Characteristics are obtained for the total area of 13.62 mln km 2 . Cluster 12 (deep oceanic areas) is the largest, 2 mln km 2 (14.7 %); it includes depressions which depths, Bouguer anomalies, and Love waves have maximum values, while the sedimentary cover thickness values are average (3000 m). Cluster 10 (transition zones and superimposed structures) has the minimum square area, 37000 km 2 (0.3 %), and maximum values of released seismic energy.
It is not feasible to obtain such zonation by a manual classification of the areas by one or several parameters. Moreover, impacts of some parameters, such as the floor relief and heat flow, are more considerable than effects from the others.
The stable parameters combinations are grouped in the XY-space to occupy vast areas, and each of such areas is almost uniform in terms of the cluster type within its limits. An exception is Cluster 7 (Fig. 13) wherein the main variations are caused by the extensive highproductivity plume magmatism (resulting in an intensive mosaic-pattern magnetic field) and the specific topography of aseismic oceanic uplifts. The described classification technique provides for differentiating between spacious structures characterized by similar parameters. A lengthy structure, such as the Gakkel ridge (the area of ultra-slow wave velocities), which width is less than 25-30 km, i.e. 15-17 km on average, may be skipped by the classification algorithm using the 25×25 km parameterization cell. In our opinion, at least two or three cells should be considered to ensure a reliable definition. Besides, an important factor is that in the conditions of slow spreading, low-productivity magmatism and poorly contrasting topography (in compa-  rison to structures of other MOR segments), structural inhomogeneities in the crust and upper mantle are less numerous than in the high-velocity conditions. The parameters and potential fields transformants, which are used in calculations, are almost identical in values at both flanks of the ridge, and the abyssal part of the Eurasian basin looks uniform in the classification. The only parameter with a contrasting value at the Gakkel ridge is the seismic moment. However, in the normalized space of parameters (see Fig. 12), no extreme seismicity values (average reference magnitude ~3) are estimated along the ridge axis, while strong seismic events are noted in the Novaya Zemlya. This could be adjustable by applying weight factors, but this kind of adjustment in the classification procedure must use justified quantitative criteria so as not to lead to purposeful fitting of the results as desired. Uniform parameter weights are used in our analysis, because the weight factor selection approach has not been justified yet. The resultant stable classification is different from any classification based only on visual analysis of the topography (the shelf areas look by far less contrasting if only their relief is considered, but its statistical correlation with other parameters reveals more contrasting patterns) -as a parameter describing the object's T a b l e 2. Square of cluster areas Т а б л и ц а 2. Площади, занимаемые кластерами  . 13. Geodynamic zoning of the Arctic region as a result of the cluster analysis, and cluster group for the main structural zones. The shelf edge is shown by the dashed line.
geometry, relief is more comprehensively represented with account of other characteristics related to energy release and the geometry of the internal boundaries of the crust and the upper mantle. The basis for the geodynamic classification is thus provided. It is also noteworthy that the visual correlation compares parameters against clearly defined extreme values of one parameter or another, while the quantitative correlation provides for comparing different background values (average values of individual areas) that are practically undetectable for the visual analysis, but very important for obtaining characteristics of spacious areas.

DESCRIPTION AND GEODYNAMIC MEANING OF THE GROUPS
Cluster analysis, as discussed above, is aimed at identifying objects that share common properties as shown by combinations of parameters. For purposes of the geodynamic analysis, interpreting of cluster groups is the task which formulation can justify the models or give grounds for objections and raise new questions. Group 1 (Clusters 3, 4, 5, 6, 9, 11, and 13)shelf and continental areas Cluster 3. Its specific features are as follows: the majority of parameters of shelf areas have average values; the isostatic anomaly, Bouguer anomaly and Love waves are the lowest; and the sedimentary cover thickness is practically maximum in the region. This cluster is revealed in the East Barents sea, Timan-Pechora and South Kara basins, and the eastern segment of the Khatanga trough (Fig. 13, Table 1). It corresponds to structures with deep positions of the sedimentary cover bottom within the shelf and continental areas. Under such conditions, propagation of surface Love waves is hindered, and the absence of positive topographic anomalies suggests a lack of mass over the compensatory isostatic surface.

6.2.1.
In Cluster 4, the majority of parameters have values typical of all the shelf areas -negative Love waves and P-waves, increased seismic moment and practically minimal magnetic field. The shelf depths are minimal in Cluster 4. This cluster as a unified whole is revealed in the area from east to west, from Alaska (Seward peninsula, and the Brooks ridge) across the Chukotka sea, the eastern Chukotka folded region, the East Siberian sea (including the Novosibirsk islands), to the edge of the Laptev sea shelf, and to the islands of Severnaya Zemlya archipelago (see Fig. 13, Table 1). It is also revealed in the southern Franz Josef Land and the northern Novaya Zemlya. Negative Love wave values may reflect the presence of the Paleozoic folded basement. Low P-wave values are due to the complex tectonic setting and high macrofracturing. Besides, this cluster comprises the zone with the quiet negative magnetic field, which is typical of the major territory of the shelf areas. Seismicity is generally higher than in other shelf zones. This cluster is consistently observed across the entire south-eastern periphery of the deep oceanic Arctic region.
Cluster 5. Its specific features are minimum values of Bouguer anomalies, Love waves and sedimentary cover thickness. The topographic and S-wave anomalies are characterized by maximum values. Other parameters have background values. In space, this cluster is adjacent to Clusters 3 and 4 (Fig. 13, Table 1) and represents a continuation of Cluster 4 in the land folded zones. This is evidenced by low values of the Bouguer anomalies, which are typical of folded structures with roots, as wells as by the reduced Love-wave velocities and the increased surface layer thickness. Positive S-wave values show that Cluster 5 is of the continental origin. The anomalous magnetic field values are above average. This cluster is revealed in the Pai-Khoi-Novaya Zemlya zone, the Taimyr peninsula, the southern periphery of the East Siberian sea, as well as in some islands of the northern Canada (which suggests the similarity of the latter areas).
Cluster 6. Its parameters are similar to those of Cluster 3 (the major sedimentary basins in the shelf areas), but the effective surface layer reconstructed from the Love-wave data is thinner, and more P-wave values are positive, which is indicative of the higher consolidation and lower tectonic fragmentation of the basement. The average sedimentary cover thickness is about 5500 m. This suggests a young age of the upper structural layer of the platform, which overlies the Paleozoic basement (Fig. 13, Table 1). This band of about 400 km is located in the northern part of the Barents sea, which Jurassic rift structures are well known, as well as in the Laptev sea and the northern periphery of Alaska and Canada.
Cluster 9 is revealed at the periphery of the Baltic, Greenland and Canadian shields (Fig. 13, Table 1). Values of the parameters are typical of the continental masses -the low Bouguer anomaly with a high average heat flow (90 mW/m 2 ), low isostatic values, rather low Love-wave values (which is indicative of a thicker effective layer), and positive S-wave values. Sufficiently high P-wave values may evidence either the mass consolidation or the general state of stresses in the crust and the lithosphere within this cluster. The sedimentary cover is 6900 m thick. The presence of hydrocarbon deposits in these areas seems highly probable, which is supported by the fact that gas and condensate fields have been discovered in the Barents sea.
Cluster 11 is revealed in the Baltic, Greenland and Canada ancient shields, the zone of the Timan-Pechora baikalides, and the Caledonian Svalbard (Fig. 13, Table  1). Its characteristics are very similar to those of Clus-ter 9, but it is distinguished by lower heat flow values, thinner sedimentary cover (to 1400 m), positive elevation values, and values of the anomalous magnetic field which are higher than typical of the shields. The spatial patterns of this cluster and other clusters of the shelf and continental areas are rather compact, without any scattering of the clouds of values, and zoning is thus quite reliable.
Cluster 13. Its characteristics are similar to those of Cluster 11, but there are significant differences (Fig. 13,  Table 1), the main of which is an increased average value of the anomalous magnetic field (310 gammas). This parameter is the main feature of Cluster 13. Other specific features are the thicker sedimentary cover, the deeper relief, and higher values of Bouguer anomalies. This cluster has a scattered spatial pattern, but is not revealed randomly everywhere. It is observed on the shields, in the White Sea, in the northern Kara Sea, and (which is very interesting) in the symmetrical frame of the Amerasian basin at the junctions of the Alpha and Mendeleev ridges and the shelf areas. This may suggest a similar origin of the magnetically active layer in the area of this cluster. It should be noted that the main parts of the basin and the Alpha ridge do not fit to this cluster.

Group 2 (Clusters 7, 12, and 14)deep oceanic areas
Cluster 12 is revealed in the major part of the Arctic Ocean, including the deep Amundsen, Nansen and Canada basins, and a significant part of the Amerasian basin (Fig. 13, Table 1). This cluster is characterized by the following maximum (or increased) values: Bouguer anomalies (246 mgal) (solid lithosphere), isostasy (45 mgal), seismic moment, relief depth (3415 m), and Love-wave velocities. Its parameters correspond to an average heat flow for the World Ocean. The sedimentary cover thickness amounts to 3,000 m. Except for the last parameter, all other parameters, including negative S-wave values, are typical of an oceanic object.
Cluster 7. Many of its parameters are similar to those of Cluster 12. However, it is distinguished, first of all, by a mosaic pattern and high values of the anomalous magnetic field in the Amerasian basin, east of the Lomonosov ridge (Fig. 13, Table 1). According to [Sokolov, 2009], this area is an extensive evidence of pluming magmatism in the absence of inversions in the magnetic field in the Cretaceous, from 120 to 80 Ma. This area is aseismic (while seismicity in the basins is medium); its average relief is slightly higher; and the average Bouguer anomaly value is decreased to 186 mgal. This combination is typical for the oceanic crust that forms due to a highly productive magmatic source [Sokolov et al., 2008]. It should be noted that Cluster 13, which magnetic field is strong, is adjacent to Cluster 7 near the continental margins. This may be indicative of a similar genesis of the crustal segments with strong magnetic anomalies -a large mantle plume may propagate underneath this territory including the shelf and continental areas.
Cluster 14 is revealed in the Norwegian-Greenland basin (the Mona and Knipovich ridges with adjacent deep troughs) and the western part of the Gakkel ridge (550 km north of the Lena trough) (see Fig. 13, Table  1). It is characterized by high and maximum average values of heat flow (129 mW/m 2 ), isostatic anomaly (intensive formation of an excess mass over the compensation surface), Love waves (minimum thickness of the effective surface layer), P-wave tomography (the state of stresses in the crust and the lithosphere), total seismic moment (intensive crustal accretion due to magma penetration along MOR and associated seismicity). The sedimentary cover is thin on average, which is normal for a young basin. Cluster 14 is characterized by parameters that are not typical of MOR -the anomalous magnetic field is low, and S-waves tomography shows positive ('cold') values. The boundaries of this cluster amount to one third of the length of the Gakkel ridge. The magnetic anomalies may be low because thick trapps were formed during opening of the basin [Sokolov, 2011], as reflected in the magnetic field (see Fig. 11), and the magma sources were depleted in Fe, which is evidenced by the low Fe concentration (7 %) in basalts sampled in the region (about 9.5 % background) [Sushchevskaya et al., 2010]. The S-wave 'cold' values are difficult to explain.
These characteristics are typical of young oceanic basins which opening commenced in the Paleocene or the Oligocene.
The profile along the Gakkel ridge (Fig. 14) reflects many geophysical parameters used in our study. A 'cold' tomographic anomaly is identified at the axis of the Gakkel ridge as a 'slab' dipping to the west from the depths about 700 km below the Laptev sea to the boundary; the 'slab's centre is underneath the central part of the ridge. It may be a relic resulting from opening of the Amerasian basin, in accordance with the rotation hypothesis [Khain, Lomize, 1995]. The western segment of the ridge is marked at the surface by a 'cold' anomaly, which is absolutely not typical of the MOR system of the Earth. This setting observed in the ultraslow spreading conditions is globally unique. There are grounds to suggest that, in the absence of pluming, spreading is not initiated by injections of hot substance, but vice versa, rifting triggers a compensating uplift of the matter and subsequent occurrence of 'hot' anomalies in the originally 'cold' blocks and facilitates spreading. In the profile, configurations of the 'slab' anomalies support the hypothesis that the slabs result from thrusting of the lithospheric blocks (according to the Arctic rotation hypothesis), rather than underthrusting or subduction. Seismicity along the Gakkel ridge and along this profile is typical of MOR with increasing depths (to 40 km) of hypocenters, while seismicity across the continental areas is scattered. In case of a weak seismic event (~3 Mb), its focal depth cannot be precisely determined, and such an event is generally assigned a numerical value that does not show a depth and refers only to the depth class marked by the same digit. For example, '10 km' refers to the shallow-foci (crustal) class, and '33 km' refers to the deep-foci (upper mantle) class. For the Gakkel ridge, values at these levels amount to almost 95 % of all the seismic events, and the remaining values refer to rare strong and deep (below 33 km) events that may occur within the MOR areas. Besides, some deep-foci events shown below the deep-ocean part of the profile may be a result of estimating the focal depths from low-precision data, which were not filtered out because actual events were shown in the eastern part of the profile. The 1999 earthquake cluster at the Gakkel ridge between 85° and 87° N is coincident with a sharp negative isostatic anomaly. This suggests rather intensive local uplifting of the matter, although an uplift and local anomalies are not contrastingly reflected in the tomography field. Central-type volcanic structures are shown in the de- Рис. 14. Профиль вдоль хребта Гаккеля с пересечением Евразии -тепловой поток, изостазия, сейсмичность, аномалия Буге, сейсмотомографический профиль по S-волнам. tailed relief scheme of the ridge in [Michael et al., 2003]. In the transition area from the tomographically anomalous 'cold' ridge segment to the 'hot' segment, the heat flow is increased to 500 mVat/m 2 and above. It is most likely that this peak can be explained by comparing the heat flow in the Spitsbergen archipelago and the anomalies in the Barents sea, rather than by analyzing the values of this parameter along the profile. Heat flow anomalies in the western Arctic (see Fig. 7) are clustered in the territory comprising the ridge's flanks and the continental margins to form a zone superimposed on the oceanic and continental structures, i.e. it does not correlate directly and everywhere with spreading. In the eastern part of the ridge, the heat flow is poorly correlated with the P-wave tomography.
Thus, Cluster 14 covers the territory with a very specific combination of parameters that are typical of the formation of a young oceanic segment in the conditions of the large transform system (represented by structures of the Knipovich ridge and the Lena trough, in our opinion) [Sokolov, 2011]. The similar scale structure is the Romanche transform fault zone in the equatorial Atlantic Ocean, but it is three times as old. This zone has been identified from the same sets of geophysical parameters in the cluster analysis study of the Atlantic Ocean, and there is a close agreement between the results in [Sokolov et al., 2008] and the cluster analysis results concerning this zone in the Arctic study. The results concerning the Greenland continental margin are also coincident.

Group 3 (Clusters 1, 2, 8, and 10) -transition zones and superimposed structures
Cluster 1 is revealed on the continental slopes of Canada and the western Barents sea, and the Podvodnikov basin between the Lomonosov and Mendeleev ridges (see Fig. 13, Table 1). In this cluster, the thick (8000 m) sedimentary cover is located on the matter characterized by 'oceanic' values of such parameters as the Bouguer anomaly, isostasy and Love waves. This cluster is distinguished as a subtype of the abovedescribed deep-ocean settings with the increased thickness of the sedimentary cover. This feature is much more significant in the Arctic than in the Atlantic Ocean. It is of interest that this cluster is present near the continental slope of the East Siberian sea in direct contact with Cluster 7, which, in our opinion, is an indicator of plume magmatism. This setting is favorable for the occurrence of hydrocarbon traps, like those known in the Norwegian sea.
Cluster 2 is distinguished by parameter values close to 'oceanic' ones, except for S-waves. However, it has an unusual combination of maximum average heat flow values (263 mW/m 2 ) and low heat flow values in the 'cold' blocks of the lithosphere, as shown by the tomog-raphy data (see Fig. 13, Table 1). It is the most interesting cluster. As mentioned above, thermal anomalies initiated by rifting are present in this area in combination with the tomographically 'cold' blocks of the lithosphere. They are located symmetrically with respect to the MOR axis and directly on the axis of the Gakkel ridge. A similar combination of parameters is established for the Spitsbergen region wherein a similar-size dome is located, according to both the heat flow measurements (the heat flow is 10 times higher than the background values) [Khutorskoy et al., 2009] and data on the Quaternary volcanic activity in the archipelago. Another thermal dome is located in the southern Barents sea. In our opinion, this setting is a good illustration of the idea that the departure of continental masses from each other can initiate compensatory ascending mantle flows that are accompanied by an increase in heat flow, and the lithosphere is thus gradually heated as shown in the S-wave images.
Cluster 8 is revealed in the continent-ocean transition zones, the Lomonosov and Mendeleyev ridges, and the Chukotka plateau (see Fig. 13, Table 1). Its main features are the 'continental' value of the Bouguer anomaly, the isostatic 'oceanic' value, and a rather large average value of the sedimentary cover thickness (3200 m). These characteristics are typical of the continental margin wherein the sedimentary bodies are growing horizontally and excessive masses are accumulated over the compensatory isostatic alignment surface. Similar parameters are typical of the Mendeleev ridge and Chukotka plateau, but absent in the Canadian coastal area. This is indicative of the presence of the continental crust behind the shelf edge in this area, at a distance up to 300 km to the deep oceanic area [Sokolov, 2009].
Cluster 10 is revealed only in areas with extreme values of the total seismic moment, specifically in the Novaya Zemlya region, the Gakkel ridge (in the area of the 1999 earthquake cluster), and the Molloy ridge (see Fig. 13, Table 1). Other parameters have average values in different zones and considered insignificant. This cluster has a scattered pattern and occupies only 0.3 % of the total square area.

CORRELATION BETWEEN THE CLUSTER ANALYSIS RESULTS AND THE BASIC GEODYNAMICS CONCEPTS FOR THE ARCTIC REGION
The established set of clusters, i.e. stable combinations of the geological and geophysical parameters, reflects the geodynamic and structural properties of the main zones in the deep Arctic and the neighboring shelf areas. It is consistent with the theoretical geodynamic history of the Arctic [Khain, 2001]. According to the well know concepts of the Arctic region evolution, opening of the deep oceanic area proceeded in two stages -the Amerasian basin opened in the period from 157 to 140 (120) Ma (or to 80 Ma), and the Eurasian basin opened 56 Ma ago. Our study provides some additional data for these concepts. In our study, the shelf areas are differentiated by combinations of the parameters, including the sedimentary cover of varying thickness, the basement with varying degrees of tectonic macrofracturing, various heat flow and seismicity, and the presence of magnetized bodies.
A more important result is variations in the most common Cluster 12 (deep oceanic areas) (see Fig. 13). The scattered Cluster 7, showing the high magnetic field superimposed on the basins, suggests that the Amerasian basin opened when pluming took place underneath the spreading axis and created local crustal structures, and the basalt layer thickness increased due to high-productivity magmatism. The symmetric relief structures [Sokolov, 2009], having numerous analogues in areas with proven presence of plumes underneath the spreading system, support the rotational hypothesis for this part of the Arctic. It is most probable that the west-dipping 'slab'-type anomalies revealed by the S-wave records underneath the Gakkel ridge axis show that the continental drifting at the Arctic periphery resulted from independent movements of the continental masses and thrusting onto the oncoming substance, and it is unlikely that the blocks were pushed apart by pluming in the basin's centre. Across the entire area, the plume markers are significantly less wide than the continental blocks involved in drifting.
Another variation in the deep-oceanic area cluster characteristics is the zone near the Mona, Knipovich and western Gakkel ridges, which is consistently identified as a particular type in the classifications for both the Arctic and the Atlantic. It is marked by the 'cold' subsurface S-wave anomaly and the positive P-wave anomaly (i.e. the state of stresses in the lithosphere), increased heat flow, and seismicity. These characteristics are consistent with the concept that the Eurasian basin opened due to spreading, and suggest that the opening resulted from own movements of the continents which facilitated spreading processes, and it is unlikely that the continents were pushed apart by ascending hot mantle flows. In the course of pluming, deep heat and magma were taken to the surface with some delay and provided for accretion of the oceanic crust. The thermal domes are observed in couples located symmetrically with respect to the MOR (which confirm the horizontal-type origin), and a few are imposed on the continents. In any case, these phenomena are not revealed along the entire opening axis and were thus unable to push apart the continental blocks. The above-described combination of the parameters is typical of the junction of the Arctic and Atlantic segments of the World Ocean along the transform fault system, which length is almost 1100 km. A similar junction of the large oceanic segments is located only in the Romanche fault zone. Discussing the origin of forces pushing apart the continental blocks is outside the scope of this paper.

CONCLUSIONS
1. The cluster analysis techniques applied to zone the crust and upper mantle of the Arctic region reveals the following stable groups supported by the geological features: • Group 1 including seven clusters (Clusters 3, 4, 5, 6, 9, 11, and 13) -shelf and continental areas; it shows significantly inhomogeneous combinations of the sedimentary cover and the basement reflected in the analyzed parameters; • Group 2 including three clusters (Clusters 7, 12, and 14) -deep oceanic areas; it shows the characteristics of abyssal areas, superimposed magmatic structures and the transition zones between the major oceanic segments • Group 3 including four clusters (Clusters 1, 2, 8, and 10) -transition zones and superimposed structures; it differentiates the ocean-continent junctions and shows specific imposed conditions in the local areas.
2. The shelf zones are differentiated by several stable combinations of parameters, including the various sedimentary cover thickness, the basement with varying degrees of tectonic macrofracturing and various states of stresses, heat flow, the presence of magnetized bodies and gravity anomalies reflecting the geometry of the lower structural layer of the crust and upper mantle.
3. In the deep oceanic areas, the cluster analysis reveals combinations of blocks with anomalously 'cold' Swave values and increased heat flow, seismicity and Pwave anomalies, which are consistent with the cluster analysis results for the Atlantic. Such a combination of parameters is revealed in the vicinity of the transform fault zone representing the junction of the Atlantic and Arctic segments of the World Ocean and suggests passive spreading that forms the new oceanic crust in the available space.
4. The superimposed thermal domes are located symmetrically with respect to the MOR axis, which shows the continental blocks divergence, but does not explain the origin of forces that lead to spreading. Thermal domes associated with this system can also occur on the adjacent continents.
5. The margins of the Eurasian continental shelf and the Lomonosov ridge are a symmetrical pair in the same cluster with continental characteristics. Similar characteristics are revealed along the transition zone to the north of the East Siberian sea, but absent along the Canada coastal area.

ACKNOWLEDGEMENTS
The study was conducted under the theme "Estimations of relationships between the floor relief of the Atlantic and western Arctic oceans, deformations of the sedimentary cover, degassing, geological hazards and the geodynamic state of the crust and upper mantle" (State Registration No. 01201459183) and supported by the Russian Foundation for Fundamental Research (Grants 15 05-05888, 13-05-12076 ofi_m, and 14-05-00122), Programs No. I.18P, I.43P, and II.3P of the Presidium, and Scientific School Program NSH_5177. 2012.5.