Joint Interpretation of Magnetotelluric and Potential Field Data From North-Eastern Norrbotten, Sweden

Potential field data in databases of the Geological Survey of Sweden (SGU) combined with newly acquired broadband magnetotelluric data are used to map and interpret geological units and structures of a 200 km by 250 km area in the Paleoproterozoic Norrbotten ore province (northern Sweden, latitudes 66°–68.5° and longitudes 19°–24°). In order to achieve this, a new approach is proposed with respect to extracting and analysing the possible correlation between modelled physical properties as well as their patterns with respect to depth variation within the crust. In this study, we propose the use of a neural net self-organising map procedure (SOM) for simplification, data reduction, and domain classification of the models derived from independent 3-D geophysical inversion of magnetotelluric, gravity, and magnetic data. The crustal model of the electrical conductivity structure was obtained from previous 3-D inversion of the magnetotelluric data. Processing and 3-D inversion of the regional magnetic and gravity field data were performed using an open-source object-oriented code called SimPEG. The input data to the SOM analysis contain resistivity, magnetic susceptibility, and density model values within the Norrbotten area for some selected depth levels of the entire crust. The domain classification is discussed with respect to the geological boundaries and composition of the crust. Consistency between model domain classification and geological boundaries is observed in general but an apparent discrepancy is noted for some areas. The reason for the apparent discrepancy is likely related to that most geological boundaries represent surface features whereas the geophysical data includes information at depth.


Introduction
An important step towards improving utilisation of geophysical methods is on integration of information obtained from several methods. Our object of interest, the Earth, is structurally complex and 3-D modelling is in most cases required. A geophysical investigation is not an end-use in itself; it rather provides a basis for subsequent interpretations and analyses. The integration of information from several methods can reduce uncertainties in model appraisal. In this paper, we show how an unsupervised classification referred to as Self Organising Map (SOM) can be used to gather information in a form that gives direct, comprehensible, and intuitive access to all generated knowledge obtained from 3-D data inversions. We illustrate the methodology by using various geophysical data from northern Sweden.
The Norrbotten region in northern Sweden as a part of the Fennoscandian (Baltic) Shield is one of the most active mining areas in Europe and hosts several of the economically important mineral deposits including precious and base metals; i.e. the Kiruna and Malmberget iron oxide-appetite (IOA) ores, the Aitik porphyry copper ores, and the Nautanen iron oxide-copper-gold (IOCG) mineral deposits. Most large mineral systems are associated with geological structures (e.g. shear zones, folds, or faults) that provide a physical pathway in order to transport metals to the surface (Huston et al., 2016). However, several factors make it difficult to identify these structures (e.g. reactivation or deformation of the rocks over time). Multiple regional geophysical surveys such as potential field and magnetotelluric (MT) data in combination with other geoscientific data can be used to provide information on the crustal and upper mantle structures. Therefore, knowledge about regional-scale structures is of importance for an optimised search for mineral deposits. MT as a deep-probing passive electromagnetic technique is used to reveal the resistivity or electrical conductivity of the subsurface that can be interpreted to represent different geological units. Graphite or carbon film, fluids, interconnected metallic minerals, or any combination of these are the most common causes for enhanced conductivity in the Earth's crust (Glover, 1996;Selway, 2014).
Several MT surveys have been widely applied from mineral exploration to lithosphere studies in Fennoscandia (e.g. De los Angeles Garcia Juanatey et al., 2013;Korja, 2007;Korja et al., 2002;Rasmussen et al., 1987). Crustal-scale modelling and interpretation of regional-scale magnetotelluric data from the Archaean and Paleoproterozoic northern part of Sweden were recently presented by Vadoodi et al. (2021) with new dense measurements which cover an area of 200 km 9 250 km. The relatively dense site spacing provides superior resolution in comparison to an earlier MT study in the area (Cherevatova et al., 2015). The interpretation focused on verification of hypothesised genetic links between occurrences of crustal-scale conductive anomalies and locations of near-surface mineral occurrences and deposits which has been noted from several MT surveys around the world (e.g. Comeau et al., 2021;Heinson et al., 2021;Sheng et al., 2021). Some of the observed middle and lower crustal high conductivity structures correlate well with the location of major shear zones, that based on geological evidence from other and similar studies elsewhere (Heinson et al., 2018;Jiang et al., 2019), may have provided the fluid path for the metal sources at many locations with mineralisation. Possible correlations of modelled electrical conductivity structures with regional magnetic and gravity anomalies were also discussed based on visual inspection of the potential field data, and the investigated area was subsequently divided into separate domains based on model and data characteristics. A similar approach was performed by Motta et al. (2019) for the Gawler Craton in Australia.
Joint interpretation based on information derived from different geophysical methods is, in general, acknowledged as important for improved and optimum geological interpretations. The basic requirement is here to ensure that the interpretations are consistent with all available geophysical data and any other information of relevance such as borehole information and surface geology. An additional requirement is to seek models that are as simple as possible at the same time as being consistent with the data. This simplicity requirement is typically implemented by using an objective function that includes a regularisation term describing the model complexity (e.g. model roughness) when applying inversion based on unstructured meshes (Constable et al., 1987) or alternatively by using a minimum number of bodies in parametric inversions (McMillan et al., 2018). The approach presented by Fullagar and Pears (2007) provides a method which combines unstructured meshes with a limited number of bodies iteratively. Phillips and Simpson (2015) provided another approach in which inversion results based on unstructured meshes are simplified by applying a model transformation referred to as 3-D terracing. Astic and Oldenburg (2019) proposed inclusion of a dynamic Gaussian mixture prior to the inversion algorithm that allows incorporating petrophysical and geological constraints.
Various approaches are proposed to tackle the integration of different geophysical methods. Joint inversion approaches that combine different types of geophysical data having different sensitivities are of particular interest. This can provide benefits or more reliable information and improve the inversion models by improved recovering of shape and properties of the structures. Some recent reviews have been published on joint inversion approaches in various fields including mineral exploration (Lelièvre & Farquharson, 2016) and lithospheric imaging (Afonso et al., 2016). In addition, the cooperative or constrained inversion approach is discussed where a single geophysical data is inverted by including information as constraints from inversion of other datasets. In general, the goal of both approaches is to improve the resolution of the final model and to provide a more consistent interpretation of Earth properties. Moorkamp (2017) provides a thorough review of different approaches and discusses and illustrates strengths and weaknesses based on a number of case studies. It is noted that incorrect prior information introduced by implementing parameter relationships could lead to biased results.
The current study complements results from a previous study (Vadoodi et al., 2021) by a combination of different geophysical data including magnetotelluric and potential field datasets. In particular, the objective is to present a method for preparing input data to the classification algorithm that allows the inclusion of depth variations of models to be taken into account. The new integration approach involves a quantitative analysis integrating the petrophysical properties associated with all three data types (magnetotelluric, magnetic, and gravity field data). The quantitative analysis is achieved by performing individual 3-D inversions of the potential field and MT data followed by unsupervised classification of petrophysical properties from the three models jointly. The unsupervised neural net self-organising map (SOM) algorithm introduced by Kohonen (1982) is utilised for the classification. Input data features to the SOM are defined such that information on depth variation of each petrophysical property enters the classification jointly with the actual values. The results of the SOM analysis are used to provide insight with respect to the correlation between the petrophysical parameters, which potentially could be utilised in subsequent joint inversion such as presented by Paasche et al. (2010) who incorporated a modified fuzzy c-means cluster analysis algorithm in cooperative inversion.
The result of the SOM analysis is used to provide a domain classification of the investigated region which then is evaluated in a regional-scale geological context.

Geological Background
The Fennoscandian (Baltic) Shield covers a large part of northern Europe and consists of the Archaean Craton in the northeast, the Paleoproterozoic Svecofennian orogen in the central part, which is bounded by the Caledonian mountains in the northwest, and the Scandinavian domain in the southwest (Sveconorwegien orogen), (Gaál & Gorbatschev, 1987) (Fig. 1). The northern part of the Fennoscandian Shield is underlain by a continental nucleus of Archaean ([ 2.5 Ga) granitic, tonalitic and amphibolitic gneisses (Gaál & Gorbatschev, 1987;Weihed et al., 2005). In northern Sweden, the exposed Archaean rocks cover only a small area within the Norrbotten Craton which was dispersed and reassembled under multiple rifting events and several episodes of subduction-related magmatism during the Paleoproterozoic . The Archaean basement is in most areas overlain by a Paleoproterozoic cover of Karelian (2.5-2.0 Ga) (related to the Paleoproterozoic rifting of the Archaean Craton) and Svecofennian (1.9 Ga) rocks (Fig. 1a). The early Svecofennian cycle of magmatism formed 1.90-1.88 Ga bimodal volcanic successions and related sedimentary rocks (Bergman et al., 2001;Weihed et al., 2005) in response to northeast-ward subduction along the rifted margin of the Archaean Craton. Two co-magmatic plutonicvolcanic rock units generated during extension-related magmatism are assigned to the Haparanda suite (HS)-Porphyrite group (PorphG) and Perthite-monzonite suites (PMS)-Kiirunavaara group (KiirG), (Bergman et al., 2001). The late Svecofennian magmatic cycle (1.81-1.62 Ga), resulted from eastwarddirected subduction with related E-W-directed crustal shortening (Bauer et al., 2018), generated widespread S-type intrusive rocks (Lina granite) associated with extensive I-to A-type felsic batholiths (marked as granites, granodiorites, syenitoids in Fig. 1a) of the Trans-Scandinavian Igneous Belt (Weihed et al., 2002).
Paleoproterozoic rocks in northern Norrbotten record at least two compressional deformation events and associated metamorphism related to accretion processes, at the early stage at 1.89 Ga and later during the 1.81-1.78 Ga period (Bauer et al., 2011;Lynch et al., 2014). The metamorphic evolution of the northern Norrbotten is poorly constrained, Vol. 179, (2022) Joint Interpretation of Magnetotelluric and Potential Field Data Raahe-Ladoga zone. The red polygon in the insert map indicates the study area however, the metamorphic grade has been suggested to increase from west to east from upper greenschist to amphibolite facies (Bergman et al., 2001).

Magnetotellurics (MT)
Broadband magnetotellurics (BMT) data were collected at a total of 104 new stations during 2015-2018 (Vadoodi et al., 2021) with 5-10 km station spacing in the frequency range of 0.001-300 Hz. Data were recorded using the MTU2000 system (Smirnov et al., 2008) and instruments manufactured by the Metronix company. Recordings were acquired with two sampling frequencies, 20 Hz for 1-day measurements while 1000 Hz was used simultaneously for 2 h during local night-time. The data processing was done by using a remote reference technique together with a robust statistical tensor estimation procedure (Smirnov, 2003). A total of 165 stations (red coloured circles in Fig. 2) including the new measurements together with previously available MT data were analysed to reveal major crustal conductors in the study area (Vadoodi et al., 2021), and the final 3-D resistivity model is used for further investigations in the current paper.

Airborne Magnetic and Ground Gravity Survey Data
The magnetic anomaly data ( Fig. 2a) selected for 3-D inversion are obtained from the 2 km by 2 km Circum-Arctic magnetic anomaly grid (CAMP-GM: Figure 2 a Map of total magnetic anomaly field from CAMP-GM: Circum-Arctic magnetic anomaly grid after 2 km upward continuation. b Bouguer anomaly data measured by SGU. The potential field data are shown as colour images and by applying shading to the data. Coloured circles represent MT stations used for inversion in the referenced study. Black-coloured symbols mark mineralisation at Aitik, Nautanen, Malmberget and Kiirunavaara deposits. Solid lines (white colour) delineate main deformation zones (DZ) in the area, (see Fig. 1) Vol. 179, (2022) Joint Interpretation of Magnetotelluric and Potential Field Data Circum Arctic Mapping Project-Gravity and Magnetic) constructed by compiling data from the world digital magnetic anomaly map (WDMAM, Hemant et al., 2007;Maus et al., 2007) and using airborne magnetic data for short-wavelength components (Gaina et al., 2011). The original magnetic data utilised in the construction of the map in Sweden are the regional-scale aeromagnetic data with survey line spacing of 200 m acquired by the Swedish Geological Survey (SGU). The total magnetic anomaly data indicate complex features such as complicated deformation zones, folded features, and other valuable information in the study area. The regional gravity measurements were made as a ground survey with an average point spacing between 1000 m and 1500 m by SGU between 1960 and1985. Gravity data are presented in the form of the Bouguer anomaly map (Fig. 2b). The high values typically coincide in general with the high-density intrusive rocks and the low values refer to areas dominated by supracrustal rocks.

MT
A three-dimensional model of electrical resistivity in the study area was provided by Vadoodi et al. (2021) by using the ModEM code (Egbert & Kelbert, 2012). For the final model, the input data consist of impedance tensor data from 165 sites in the period range of 0.003-1000 s (4 periods per decade) and tipper data from 27 sites (in period 1-1000 s, 4 periods per decade). The model domain was discretised into 153 9 126 9 54 cells with horizontal cell size 2 km 9 2 km. Padding cells increased gradually in size by a multiplication factor of 1.5. The uppermost cell thickness was chosen as 50 m and the cell thickness increased downwards by a multiplication factor of 1.2 to a total depth of about 5700 km. The model domain covers 8900 km in the north direction and 8800 km in the east direction. Several choices of regularisation parameters, data partitioning and starting models were tested prior to defining parameters for the final inversion. The data were fitted to a total RMS (root mean square) value of 2.1 when using a 5% error floor for the impedance tensor elements and an absolute error of 0.05 for tipper data. For more details on inversion strategy and data misfits, please see Sect. 4 in Vadoodi et al. (2021). Figure 3 shows the final 3-D resistivity model for six depth levels selected for further analysis jointly with density and magnetic susceptibility models described in the following Sect. 5. The 3-D resistivity model indicates a highly heterogeneous distribution of conductivities with very high conductance of more than 3000 S at some places in a generally resistive crust.

Potential Field
The potential field inversion models were obtained using an open-source object-oriented code written in Python called SimpPEG (Simulation and Parameter Estimation in Geophysics, Cockett et al., 2015). SimPEG includes various facilities for inversion such as staggered grids, standard numerical solver packages, optimisation algorithms, and several choices with respect to model parameterisations, a priori models, and smoothness constraints. Since we focus on regional structures in the study area, an upward continuation filtering of 2 km for magnetic data and 500 m for gravity data was performed before inversion.
The magnetic data inversion performed for this study assumes that magnetisation is induced (no remanent magnetisation) and that the main inducing magnetic field has a field strength of 51,700 nT, a declination of 3.7°, and an inclination of 76.3°t hroughout the model domain (International Geomagnetic Reference Field, IAGA Working Group 1, 1981). An octree-based discretisation is defined for the model. The core region covered by the data is discretised with a 1.5 km horizontal resolution and a 300 m vertical resolution. The choice of model resolution is mainly governed by the use of upward c Figure 3 Sections from 3-D resistivity model (Vadoodi et al., 2021) at different depth levels. a 5 km, b 10 km, c 15 km, d 20 km, e 25 km, f 30 km. MT site locations are marked with white-colour circular symbols continued data in the inversion but is also chosen to be compatible with the discretisation used for the MT model. Padding distance around the core region is defined to 15 km in north and east directions. The model domain covers 380 km in both north and east directions and 120 km in vertical direction. Topography varies gently in the range of 100-500 m. Topography was assumed flat in the inversion since the topography variation in the overlapping MT model area is limited to ca. 200 m and thereby insignificant topographic contributions at the 2 km upward continuation height. Model susceptibility was constrained to be in the range of 0.00001-2 SI units and a background and starting susceptibility of 0.001 SI was used. The data uncertainty is set to 3 nT for all magnetic data points. The total number of magnetic input data was 2891 with a regular distribution and data point spacing of 5.1 km. The model discretisation for Bouguer gravity data is defined similarly to that for magnetic data based on an octree mesh. The core cell size is defined as 1 km in horizontal and 500 m in vertical direction. Padding distance around the core region is defined as 10 km in horizontal directions. The dimension of this model domain is 512 km in north, 510 km in east, and 120 km in vertical direction. Data uncertainty is assumed to be 0.75 mGal. The total number of gravity input data was 3751 with a regular distribution of 4 km data point spacing. A Tikhonov regularisation was used in the inversion to impose smoothness constraints for both gravity and magnetic data inversions. Inversion is performed by a Gauss-Newton (GN) approach with an optimisation using the Conjugate Gradient (CG) method. A maximum of 30 GN iterations and a maximum of 10 CG iterations with CG tolerance 0.001 were set. Figures 4 and 5 present horizontal depth slices for the final density and magnetic susceptibility models, respectively. A RMS value of 0.8 mGal and 3 nT was achieved for the gravity and the magnetic data inversion, respectively. Figure 6 presents the normalised data residuals for magnetic and gravity data points after interpolation to grids. Calculations are based on residuals between the observed and predicted data, normalised by the data uncertainties. The grids are shown using a linear colour scale. Calculated standard deviation was 0.88918 and median 0.030 for gravity data. Standard deviation was 0.92130 and median 0.049 for magnetic data. The lack of gravity data in the eastern part of the area resulted in an incomplete grid.

Self-Organising Maps (SOM)
The self-organising map algorithm developed by Kohonen (1982) is an unsupervised clustering method using a type of neural network algorithm to simplify analysis, interpretation, and visualisation of a high-dimensional dataset. In short, the high-dimensional data are projected in a non-linear manner onto a two-dimensional grid composed of hexagonal cells (nodes) referred to as the SOM space, such that similar input data appear adjacent to each other. Below, the principles of the SOM are briefly introduced. For more information, the reader is referred to Kohonen (1982), Vesanto et al. (2000), and Fraser and Dickson (2007).
In a SOM analysis, multi-dimensional data after normalisation of each data type are presented as n-dimensional (n-D) vectors (where n is equal to the number of input variables) and then provided as input to the SOM algorithm. The n-D vectors are seeded (typically randomly) by a defined smaller number of seed-vectors equal to the size of the required output map in the SOM space (e.g. 10 9 15 sized map means 150 seed-vectors). In an iterative two-step learning process, all seed vectors are trained to represent the patterns within the input data. First, each input sample, irrespectively of their geographic location, is compared with all seed-vectors within a particular radius of the input sample based on a measure of vector similarity, typically an Euclidean distance measure. Similar data vectors after clustering are approximated by a single vector referred to as Best Matching Unit c Figure 4 Sections from 3-D density contrast model at different depth levels. a 5 km, b 10 km, c 15 km, d 20 km, e 25 km, f 30 km. MT site locations are marked with white-colour circular symbols. Den-C density contrast (BMU). The BMUs are typically organised in a twodimensional map referred to as the SOM space such that similar BMU vectors appear adjacent to each other whereas dissimilar BMUs are well separated. The process is repeated several times for each input data vector in the data space to achieve the best approximation of a set of data vectors. Thus, after training, each BMU vector represents a group or cluster of the data vectors. The similarity between BMUs in adjacent cells is visualised in a unified distance matrix called the U-matrix plot, where colours are used to visualise distance relations between the adjacent BMUs. Simplification or data reduction of the input data is achieved since each BMU typically approximates several and similar input vectors. The geographic information of each input vector is tagged to the data approximated by each BMU and this enables visualisation of the clustering in geographic space after performing the SOM analysis. Component plots are another visualisation of the nodes representing the variation of a particular variable on the 2-D SOM space. Comparison of the values at the same location in SOM space can be achieved by visual inspection of the components map which allows detecting the relationship and correlation among all input data types.
In the current study, the SOM analysis is performed by using the SiroSOM software distributed by the Commonwealth Scientific and Industrial Research Organisation (CSIRO), (Fraser & Dickson, 2007). Since the input data contain values with different units, the values are typically normalised either by the standard deviation after subtraction of the mean or by range to obtain values between 0 and 1.
The input data vectors contain final 3-D resistivity, magnetic susceptibility, and density model parameter values presented in Sect. 4 for six different depth levels (5, 10, 15, 20, 25, and 30 km) sampled with a regular 2 km grid distance horizontally. The Normalised data residuals maps (residuals between observed and predicted data normalised by data uncertainties) from a magnetic and b gravity data. The gap on the right side of the gravity residuals map is caused by lack of data Vol. 179, (2022) Joint Interpretation of Magnetotelluric and Potential Field Data same sampling is needed for each model, so we interpolated model values by using a linear interpolation method in order to obtain coincident sample locations for the three models. Since models were obtained with inclusion of smoothness constraints in the inversion, errors introduced by interpolation are negligible in relation to the present analysis with focus on regional structures. The choice of treating data from different depth levels separately is done to include depth variation of the petrophysical model parameters in the analysis and classification. A total of 11,374 data vectors of dimension 18 (3 petrophysical properties and six depth levels) corresponding to each easting and northing location were input to the SOM analysis. The mean value and standard deviation for each data type were used for normalisation of data prior to analysis. This normalisation is preferred compared to normalisation based on data in each depth level separately since the information of the relative data variation with depth is maintained and thereby enters the analysis as data features in addition to the actual data values.
The 2-D SOM grid/matrix consists of 40 9 34 cells which corresponds roughly to a data reduction by a factor of 150: (11,374 9 18) input data/ (40 9 34) SOM output. Component maps referring to each normalised input parameter to the SOM analysis are shown in Fig. 7 together with the U-matrix. One of the cells in the SOM grid (Umatrix) is marked as a reference for illustrating the use of the component maps and how to link the data in the SOM space to positions in geographic space. The approximation of the input data by this specific BMU can be seen directly by inspecting the same location in each component plot. Thus, the resistivity at 5 km and 10 km depth are noted to be of intermediate values whereas the resistivities at the four other depth sections are low. Susceptibilities are seen to be low at all depth levels. The BMU in this example is associated with 12 geographical locations as illustrated and discussed in Sect. 5.2. The distribution function for the magnetic susceptibilities is highly skewed towards zero values and the minimum values for the legends are therefore written as zero due to truncation of the decimals.  Component maps of the SOM analysis corresponding to each normalised input parameter in the SOM space. U-matrix represents the unitless Euclidean distance to the neighbouring cells displayed in the SOM. The initial units for resistivity, density contrast, and magnetic susceptibility are Xm, g/cm 3 , and SI Units, respectively. The black circle on the top right side of the U-matrix map shows the selected cell (BMU) used to exemplify the link between data in the SOM space and geographics space 1080 R. Vadoodi and T. M. Rasmussen Pure Appl. Geophys. maps are particularly useful for visualising correlation between input data. It is evident from Fig. 7 that no simple linear correlation exists between the distribution of the three petrophysical properties, but some depth-wise correlation is noted with respect to depth variations for a single petrophysical property. The latter is expected due to the smoothness constraints in the applied inversion algorithms.

Spatial Analysis of the SOM
In Fig. 8, the classified data in the SOM multivariate data space are plotted back into the geographical space and visualised by using a unique colour coding of the SOM-matrix (Fig. 8a). Identical colours in Fig. 8b refer to data associated with the same BMU and thereby approximately similar values in all three petrophysical properties as well as similar depth-wise distributions. Some BMUs may not be associated with data which typically occurs in the cases with distinct clusters. Linking data associated with a specific BMU to locations in geographic space is exemplified by the BMU marked with a black circle in Figs. 7 and 8a. The corresponding 12 locations in geographic space are marked with black circles on the map shown in Fig. 8b. The colour coding of the U-matrix with smooth colour variations implies that distinct contrasts of the U-matrix when Figure 8 a Unique colour coding defined in the SOM space. b The domain classification visualised by colour coding defined in a and projected back to the geographical space. White-coloured circles mark mineralisation at Aitik and Nautanen deposits. The solid lines (white colour) delineate modified main deformation zones (DZ) in the area (see Fig. 1). The thick grey line is the isoline for Nd ratio of -3 and the thin grey line isoline for Nd ¼ 0 of the LJZ interpreted as the Archaean-Proterozoic Boundary (Ö hlander et al., 1993). The black circle in a shows the marked BMU in U-matrix (Fig. 7) as an example which is associated with 12 geographical locations marked by black circles in b Vol. 179, (2022) Joint Interpretation of Magnetotelluric and Potential Field Data projected to geographic space and visualised by the colour coding refer to large differences in BMU vectors (Brethes et al., 2016). Thus, boundaries defined by distinct colour contrast can be used to define geological domain boundaries. A complementary and quantitative outlining of domain boundaries can be performed by analysing BMU-vector similarities in geographic space. The Euclidean distance between two adjacent BMUs in geographic space is calculated for each direction in the grid (easting and northing directions) and the calculated value is associated with a geographic coordinate defined by the mid-point between the two BMUs. These BMU distances from both directions are merged into a common grid and displayed in Fig. 9. Figure 10 shows the SOM colour coding superimposed on a 3-D relief map of the Euclidean distance grid. The information provided by the colour coding and distance measure is complementary. Whereas the distance measure highlights domain boundaries by showing relative differences, the colour coding when linked to the component maps and colour coding of the U-matrix provides information about absolute but normalised values of the parameters. One advantage of SOM compared to e.g. application of hard classifiers such as k-mean clustering is that the classification is able to handle and visualise data that do not necessarily have well separated and distinct clusters. The use of smoothness constraints in the 3-D inversion implies that spatial model parameter variations often are smooth and thereby violates the basic assumption entering applications involving hard classifiers. The BMU distances in Fig. 9 emphasise the edges which separate the domains visualised with different colours referring to properties within each domain. The BMUs visualised in Figs. 8 and 10 bounded by ridges defined by the distance measure have in some cases smooth variation in properties and thereby also smooth variation in colour coding within a specific domain. The latter is linked to that the BMU in the SOM-matrix are organised so that adjacent BMUs are numerically close.
The north-eastern and south-western parts of the map in Fig. 8 have larger domains compared to the other areas and this is likely a result of a smaller density of MT-stations and thereby very smooth electrical resistivity variations in the model. A good correlation between domain boundaries and some of the major deformation zones shown in Fig. 8 is noted. This is particularly the case for the NDZ where distinct differences in crustal properties are noted across the zone. In order to evaluate the applicability of the proposed analysis method, the two areas marked A and B in Fig. 9 are selected for validation purposes. Some caution must however be exercised with respect to correlation evaluation since the construction of the geological map is partly guided by the interpretation of the potential field data. The exposed areas cover only about 1% of the investigated area and the divisions on the geological map are in many cases based on inspection of the aeromagnetic data. Figure 11 shows the SOM classification and the geological map (Bergman et al., 2001) from an area hosting a number of important mineral deposits (Aitik, Malmberget, Nautanen). The area is marked A in Fig. 9. A number of major deformation zones, including the NDZ, are shown on the map as thin black lines. The SOM classification includes information about the entire crustal section and a correlation to geological boundaries based mainly on surface observations is therefore not necessarily expected. Furthermore, the geophysical models entering the classification are showing regional features only. Nevertheless, some SOM domains show good spatial correlation with mapped geological units and location of mapped deformation zones. The most pronounced discrepancy between mapped geology and the domain boundaries is noted for the large areas dominated by the Lina granites (see Fig. 1 for geological classification). A north-south trending boundary is here noted in the SOM results (marked by thick black line) entirely within the domain mapped as Lina granites. The middle crustal section towards east of the boundary is characterised by distinct lower resistivities and lower susceptibilities compared to the western side, a difference that contributes to separation into two domains. N-S elongated occurrences of meta-basalt and meta andesites coincide spatially with the SOM domain boundary, but an interpretation of any major boundary zone linked to these occurrences is not evident from prior geological studies. Figure 12 shows the SOM classification and the geological map (Bergman et al., 2001) from the area marked B in Fig. 9. The SOM classification discriminates between a number of domains and a good correlation between domain boundaries (marked with dashed black lines) and the boundaries between geological units are noted. Apparent discrepancies are however noted as illustrated by the domain boundary marked by the black line. The marked boundary can here be attributed mainly to pronounced density anomalies trending roughly N-S (see Fig. 4). The presence of SW-NE trending boundaries is linked to modelled conductive structures at middle and lower crustal depth (see Fig. 3).

Discussion
The current study presents a new way to perform joint analysis and visualisation of 3-D models with different physical properties using an unsupervised SOM classification approach. The methodology is sufficiently versatile to include more than three properties as those used in this example. An advantage of the methodology is the ability to visualise several 3-D models jointly by projection onto twodimensional maps and perform back-projection to Figure 12 The surface geology map (a) and SOM map (b) for the selected area marked B in Fig. 9. (See also Fig. 1). The solid lines show the lineaments on the surface. The bold black lines indicate the locations used for comparison with respect to the geology map geographical space. The present study describes how to include spatial variability such as variation with depth in the classification. Inclusion of variations horizontally in the classification is possible by adding e.g. gradients of model parameters as input data. The models from inversions are non-unique with respect to model properties and some care needs to be exercised in the evaluation. Nevertheless, the approach is considered valuable by providing a data reduction (model parameters) that highlights the main features of the models. No attempt was made in this study to reduce model uncertainties by adding a priori information such as petrophysical information in the inversion as described by Astic and Oldenburg, (2019), assuming relationships between the three physical properties (Lelièvre et al., 2012), application of cross-gradient techniques (Fregoso & Gallardo, 2009;Gallardo & Meju, 2003) or jointly inverting data from two methods responding to the same physical property (Blatter et al., 2019). Moorkamp (2021) presented joint inversion results of gravity and magnetotelluric data in which variation of information between density and electrical resistivity was minimised in a joint inversion based on a cross-gradient inversion scheme. The SOM analysis can obviously also be applied to models obtained from application of joint inversions. Moorkamp (2021) visualised the relationship between obtained densities and resistivities using a cross-plot (scatter-plot), which works well when only two parameters are used. The SOM approach particularly facilitates investigations of correlation between petrophysical properties when multiple parameters are included as done in the current paper, and it also provides a method to classify spatial variations. The SOM results derived from individual and independent inverted data sets may provide a first unbiased estimate of the petrophysical property distributions and provide a basis for more advanced joint inversion schemes which involves updating prior information as described in Astic et al. (2021). A potential risk with such a procedure is that information in data is used twice which in general should be avoided. The current study involves 3-D models covering a 200 km 9 250 km area and we do therefore expect that correlation between petrophysical properties vary considerably within the model space. Such a variability, which furthermore typically are very uncertain, may be difficult to build into inversion schemes based on quantification of direct coupling between properties (Moorkamp et al., 2011).
The SOM performs an unsupervised domain classification which shows areas with similar properties and similar depth variations. The classification may e.g. be used if a specific geological phenomenon such as an observed presence of mineralisation in one domain is used to predict the favourability of similar mineralisation in areas belonging to the same domain class (Stensgaard et al., 2006) provided that the domains refer to the same geological environment.
The presented visual comparison between SOM domains and mapped geological units provides examples where a clear consistency is observed but other examples are evident where more work is needed to analyse the reason for an apparent discrepancy. These discrepancies may be due to the fact that the SOM analysis includes information about depth variations whereas the geological map mainly is based on surface observations. Furthermore, the degree of exposure is fairly limited within large parts of the study area. The exposed areas cover only about 1% of the investigated area and the divisions on the geological map are to some extent guided by visual inspection of the aeromagnetic data. Validation of the methodology is therefore not straightforward. The present study is performed with a focus on mapping and interpretation of regional-scale structures and model discretisations were rather coarse. The model parameters are therefore only representative of large volumes averages and furthermore without uncertainty estimates. The result of the SOM classification should ideally be used to guide follow-up work on improved geological mapping and interpretation.
The input to the SOM analysis consists of data vectors defined at each node in the surface grid. Although the inclusion of other data such as geochemical data and petrophysical data is straightforward by expanding the dimensionality of the input vector, it is important to consider that these data types may not have any relationship to model parameters associated with the deeper part of the models. Such a lack of relationships may impact classification performance in a negative manner. Vol. 179, (2022) Joint Interpretation of Magnetotelluric and Potential Field Data

Conclusion
The SOM analysis of three petrophysical properties obtained from 3-D models has been shown to be an easy way to visualise and characterise the three geophysical models jointly to identify locations that are similar with respect to all three properties as well as similar depth-wise distribution. The component plots provide an overview of relationships between input parameters in an organised manner which combined with the U-matrix can be used to identify clustering in the data. The methodology proposed is not restricted to only three types of models which further emphasise the utility of the visualisation technique. An advantage of the SOM analysis compared to for example k-mean clustering is that it is able to characterise models having smooth spatial variations of their properties.
The clusters projected back to geographical space provide information that is directly useful for mapping and interpretation of the geology. In particular, simplifications are provided with respect to visualisation since several models are described in a single map instead of several individual maps. The results provide a basis for subsequent follow-up with respect to geological fieldwork for an integrated interpretation and analysis. The domain classification can be utilised for example in mineral prospectivity mapping or other prediction tasks where information from several geophysical methods typically needs to be integrated.
The present study included clustering with respect to depth variation. Several other spatial feature or pattern definitions in 3-D are possible and can be included in further work based on the SOM analysis (Ortiz et al., 2014).