Abstract
This study’s main objective is to better define and understand results for the most commonly used inversion algorithms in magnetotelluric data interpretation as part of geological exploration of the region of the Dolsk fault and the Odra fault. The data obtained from the eastern part of Fore-Sudetic Monocline measurements were used to describe the boundaries of lithospheric blocks (terranes) and recognize their origin. The magnetotelluric (MT) soundings were carried out to achieve this goal. There were conducted 51 soundings on five quasiparallel profiles. That allows constructing a quasiregular mesh in the area of the Fore-Sudetic Monocline. This arrangement of the measuring grid allowed reducing the influence of the largest sources of disturbances on MT data. 1D and 2D models were created by using the inverse algorithms. The models were prepared for each profile separately. Further, parallel (ModEM) 3D inversion codes were applied. The area where the investigation was done involves the region of the Dolsk fault and the Odra fault. These zones are essential geologic borders of a regional nature, and they pull apart the crust blocks with different origins. It was vitally needed to correctly identify the crust and upper mantle structure around a part of the Fore-Sudetic Monocline. The paper shows how these key features of the geological structures are revealed using 1D, 2D, and 3D algorithms.
1. Introduction
These studies’ main objective was to investigate the significance of the dimensionality of the inversion algorithm used to understand better the genesis of lithosphere blocks in the East European platforms. The data analyzed in this paper are results of magnetotelluric measurements carried out between 2016 and 2019. The measurements were made with the Geomag apparatus of Ukrainian origin by the Institute of Geophysics, Polish Academy of Sciences. Data from 51 soundings in the area of the Fore-Sudetic Monocline in Poland were analyzed and interpreted. The study area includes two significant fault zones: the Dolsk fault and the Odra fault. These are essential geological boundaries of regional character. The project’s scientific aim is to obtain a detailed image of the Earth’s crust and upper mantle structure around the northwestern part of the Fore-Sudetic Monocline. This area is a part of the Trans-European Suture Zone (TESZ), and it represents a broad and complex zone of terrane accretion at the junction of the Proterozoic lithosphere of the East European Craton (EEC) and the younger Palaeozoic lithosphere of Western and Central Europe [1, 2]. The geophysical surveys conducted in this area were mainly based on passive gravity [3], magnetometric methods, and active seismic methods [4]. These studies provided new information on the deep tectonic setting of this part of Poland [5].
Special attention was drawn to an area covering the Dolsk and Odra fault zones. Both of them are essential geological boundaries of a regional dimension. In these zones’ fields, to the southeast of the Variscan front, there is the Wolsztyn-Leszno High (WLH). The origin of the crust of this area is a disputable matter. It is located in the southwestern part of the foreland basin. A significant issue is the nature and significance of the Dolsk fault and the extent of the southwest edge of the Baltica [6, 7].
About fifty soundings on five quasiparallel profiles were conducted. All the data were collected by the Geomag Fluxgate Magnetometer apparatus. That allows constructing a regular mesh in the area of the Fore-Sudetic Monocline. Processing and preliminary data interpretation were made according to the progress of fieldwork. The analysis of dimensionality (skew, effective impedance, and horizontal magnetic tensor) shows that the studied geological structure has a 3D character, while in some regions, a situation close to 2D is observed. To compare the effectiveness, 1D and 2D models were created using the most common inverse algorithms (Occam’s and Nonlinear Conjugate Gradient methods). The models were prepared for each profile separately. Finally, parallel (ModEM) 3D inversion codes were applied. The ModEM is an inversion code that employs MPI and which, besides impedances, includes tippers. Our study’s results shed additional light on the geology and geotectonic problems of this area. Besides, an analysis of selected invariants of the transfer function has been implemented. We mainly compare the most frequently used inversion algorithms. The obtained models of conductivity distribution document deep conductive faults and show the differences in distinctive crust blocks’ resistivity structure. It allows a relevant supplement of preceding knowledge to obtain delighted recognition of the Dolsk and Odra fault zones. This research project substantially permits for better reconstruction of the tectonic evolution of the Baltica foredeep. The results also allowed assessing the usefulness of algorithms of different dimensions in the case of modeling such a complicated structure.
2. Geological Background
As it was mentioned before, the studied area is located for the most part on the Fore-Sudetic Monocline, according to regionalization from 2011 [8]. A wide terrane accretion zone separates the Proterozoic lithosphere of the Baltic Shield and the East European Craton from Western and Central Europe’s younger Palaeozoic lithosphere [2]. Southwestern Poland is a broad zone of lithospheric blocks. The northeast borders with the Proterozoic lithosphere of the Baltic Shield and the East European Craton (trans-European seam, TESZ). In turn, it borders with the Fore-Sudetic block [9].
This area is characterized by a highly complicated geological structure and a vague history of tectonic evolution. The Precambrian platform pedestal and its Lower Palaeozoic cover retreat southwest [10, 11]. In the broader sense of European geology, this means that the hypothetical Caledonian tectonic suture, which marks the collision site of Avalon and Baltica, is not located along the TTZ, but further southwest, i.e., in northern Germany and southwestern Poland [12], which is the subject of the conducted research.
According to the latest hypotheses, under the TTZ, there may be an extended and shaded border of the Palaeocontinent—Baltica [12]. That means that the TTZ formed in the Precambrian, most likely in the early Palaeozoic; it was no longer the Baltica border [8]. Unfortunately, the range of the currently mentioned Baltica rim is not fully understood. However, it can be assumed that it extends from the Dolsk fault or even the Odra fault zone [13]. Figure 1 shows the Polish Basin at the background of the main basement units, emphasizing the Variscides area, which coincides with the research area presented in this paper. The Odra fault zone is a significant discontinuity zone that separates the Lower Silesian Block from the Fore-Sudetic Monocline. It stretches for about 300 kilometers in the WNW-ESE direction, crossing the Mesozoic, Palaeozoic, and Proterozoic formations. Gravimetric and magnetic anomalies strongly mark this fault [14, 15]. The fault zone consists mainly of igneous and metamorphic rocks (both high and low metamorphism) [16].
The Fore-Sudetic Monocline substrate is a shell block located lower concerning the entire Lower Silesian Block. It was created by Permo-Mesozoic succession, which is a platform cover of the block discussed earlier [10]. It falls monoclinally at an angle of several degrees towards the northeast. The Permo-Mesozoic level rocks lie inconsistently on the folded Palaeozoic formations, which are seen in the cross-section above (Figure 2) [14].
Previous geophysical surveys in this area were mainly passive gravimetric and magnetometric methods and active seismic methods. They provided new information on the deep tectonic assumptions of this part of Poland [11, 17].
3. Method
The magnetotelluric (MT) method uses, as a source, the natural variations of the Earth’s magnetic field due to the interaction of the solar wind with the magnetosphere and ionosphere and thunderstorm activity in equatorial regions.
It is based on the electromagnetic induction process [18]. Alternating currents are induced in the rocks constituting the Earth’s interior; their directions and intensities depend on the electrical conductivity distribution [19]. Analyzing relations between the values of respective electromagnetic field components recorded at the Earth’s surface lets us conclude the basement’s electrical conductivity values [20]. It is a significant physical parameter, owing to which we can identify geological structures that differ in the electrical properties of rocks. The electrical resistivity of rocks changes in an extensive range, from a few and more (sedimentary rocks) to as many as hundreds of thousands (crystalline rocks) of ohmmeters [21]. Geophysical identification of geological structures is the basis for sophisticated geological-geophysical interpretation, assigning the separated forms to geological units and verifying the dynamical evolution scenarios [22, 23]. Electromagnetic (EM) sounding methods measure these currents directly or their respective secondary magnetic field, thus providing information on the crust’s electrical state or even the upper mantle [24]. Moreover, monitoring these parameters provides us with information about anomalous magnetic or telluric signals and changes in the conductivity distribution in the lithosphere. [25].
Magnetotelluric data processing is aimed at estimating transfer functions between the measured field components, which carry information on the size and distribution of electrical conductivity in the ground due to electromagnetic induction. Here, we are interested in the off-diagonal (, ) and the diagonal (, ) elements of the impedance tensor; also, the geomagnetic transfer function or tipper (, ) is computed from and which are horizontal components of magnetic induction; is the vertical component of the magnetic field; and and are the horizontal components of electric field which are all measured in field. These fulfil the following equation:
They are complex functions of frequency or period . The impedance elements give us apparent resistivities and and phases and . The geomagnetic transfer function is a real and imaginary induction arrow consisting of the components and each. We use Wiese’s convention [26], where real arrows point away from good conductors [27].
4. Fieldwork and Processing
Our work’s first step was to complete our data set of the study area with additional measurements. Collected data were subject to numerical processing. The most crucial step involves executing the array of about 51 long-period magnetotelluric soundings in points distributed as steadily as possible to determine the 3D model. A magnetotelluric station comprehends measurements of two telluric components in the north and east directions (electric field Ex and Ey) and three orthogonal components of the magnetic field in the north, east, and vertical downward directions (magnetic field Hx, Hy, and Hz) [28]. All measurements were carried out in a geographic coordinate system. The -axis was oriented in the magnetic N-S direction pointing to the north, and the -axis was oriented in the E-W direction pointing to the east.
All soundings were measured by the Institute of Geophysics, Polish Academy of Sciences between 2016 and 2019. These five profiles consist of typical long-period magnetotelluric deep soundings. The frequency band for this data extends from a few Hz to 0.0001 Hz. The instrumentation consisted of GEOMAG stations (based on fluxgate magnetometers of Ukrainian construction [29]). The sampling interval was 10 Hz for all stations, and they were synchronized by GPS or DCF signal. The measurement time was amounted to about ten days per station. When the recording was interrupted, the measurement was repeated for various reasons (apparatus failure, battery discharge, or intersected cables by the nearby fauna). Several soundings were duplicated due to the numerous disturbances seen in the time series [30, 31]. The studied area is challenging for the practical application of the magnetotelluric method. The main difficulties during the works included electromagnetic interference associated with numerous high-voltage power lines (in particular in the Legnica-Głogów Copper District) and interference related to DC-electrified railway lines (Figure 3.) and electrical installations of multiple localities in this area. Due to the proximity of electromagnetic disturbances, it was challenging to obtain high-quality measurement data. Polish railways (if electrified) run with direct current (DC); they cause leakage of correlated electromagnetic noise over dozens of kilometers from the tracks, which can heavily distort transfer functions [32].
The improvement of data quality required several additional measurements and calculation works. Magnetotelluric measurements were carried out simultaneously in the field and the so-called magnetic remote reference [33]. This procedure significantly reduces the interference, ultimately translating into a better look at MT curves [34].
The reference point during fieldwork was located in two places at the same time. The main point acting as a magnetic reference was a point located near the town of Birzai in Lithuania (56°1014.81 N, 24°5124.52 E), and the auxiliary point located near Suwałki in Poland (54°0044.3 N, 23°1050.8 E). We decided to conduct two soundings fulfilling the function of reference measurement, due to a long time of recording at measuring points located in the study area. Therefore, it was essential to have the same registration from the uninterrupted reference point and continuous registration at the measuring point.
5. Inversion
All data were processed using the robust remote reference method [33, 34]. The data interpretation was performed using one- (1D), two- (2D), or three-dimensional (3D) algorithms based on smooth resistivity models. The necessary foundation of 1D Occam’s algorithm is to aspire to achieve the smoothest solution [35]. To depreciate the process in 2D inversion, the NLCG (Nonlinear Conjugate Gradient) is used. This method involves the iterative matching of the specified medium to the measurement curves [36]. It leads to the gradual modification of the starting model to achieve optimal results. The WinGLink Software has been used to carry out one-dimensional Occam’s inversion and two-dimensional NLCG inversion. All data were processed by the robust remote reference method. The 3D inversion was performed based on a smooth model (ModEM). An example of individual sounding results is shown in Figure 2.
An essential element of this stage is the rotation of the soundings by an appropriate angle. This angle was determined by analyzing the polar diagrams (), phase tensor [37], and Swift skew [38] (Figure 4). Skew is one of the dimensionality data MT analyses [20]. Swift skew is a ratio of diagonal components ( and ) to the off-diagonal components ( and ). If the value of , it has a 3D character. If , the data has 1D or 2D character [39].
The analysis of polar diagrams revealed that for the central part of the research area, the structure is two-dimensional with a 42 deg strike, but three-dimensional disturbances occurred in greater depths. However, a fixed angle of rotation is assumed for the whole profile. All soundings were rotated 42 degrees so that the -axis ran along the plane of the studied fault zones.
Quantitative interpretation of MT data determines the distribution of resistivity in the geological medium and the relationship between geoelectric complexes, separated based on the analysis of measurement curves, with geological complexes. In the first step, quantitative interpretation was based on the 1D inversion method. The 1D inversion was carried out for TE, TM, and Effective Impedance (INV). All of the results were compared and contrasted, but finally, the last one is presented in this paper. It combines sensitivity to local changes of TE and regional resolution of TM. The RMS error for every sounding was reduced below 2% for the 1D inversion. In WinGLink software, the maximum number of layers was set at 45 for each sounding and was decreased only in extreme cases in a couple of soundings, but never less than 30 layers, with starting resistivity at the level of 100 Ωm. Phase and amplitude curves for one polarization (TE, TM, or INV) were inverted simultaneously. The maximum number of iterations was defined as 30 repetitions for each sounding. In one-dimensional modelling, the medium’s geometry is reduced to a flat-parallel layer system, where the resistivity only changes with depth. Inversion consists of multiple executions of a simple task, with new parameters introduced in subsequent iterations in interpretative models [28].
The next step was to perform the 2D inversion. A two-dimensional inversion algorithm consists of the iterative matching of a two-dimensional geoelectric medium to amplitude and phase measurement curves, minimizing nonlinear conjugate gradients. The error function describing the cumulative discrepancy between the empirical and theoretically calculated data for the model is reduced in subsequent iterative steps, leading to a gradual modification of the starting model. A 2D inversion was performed for several dozen initial models with variable damping parameters for the optimal solution. The answer to the inversion issue is ambiguous. Many different models can be adjusted to empirical data with satisfactory accuracy. We are talking here about the phenomenon of equivalence [18]. For the correct use of the algorithm, it is crucial to impose the appropriate mesh to make the obtained solution realistic and the proper initial conditions in the form of a starting model [31]. Such a model can be a homogeneous half-space. It is the simplest model but often gives the best results; i.e., it does not introduce artificial artefacts. For the 2D inversion, a homogeneous half-space with a resistivity of 100 Ωm has been taken as the initial model. A static shift correction has been assumed. The initial parameters were considered: a smoothing coefficient τ at 3, which seems typical for most MT inversions, while the maximum error floor was set at 5%. The number of iterations for each inversion was 100 repetitions. If the final model was not satisfactory (too high RMS or geological incompatibility), the operation was repeated, with the last model being treated as the starting model of the subsequent inversion, but with other values of smoothing coefficients (often lower than 3) and higher error floor values, but not higher than 7.5%. Due to a very rare grid, additional smoothing parameters (horizontal α and vertical β) were often used to change their values in such a way as to avoid the formation of artefacts associated with long distances between stations. The values of these parameters have not increased above 2.5. These parameters vary from about 2.5 for the rough-sampled profiles to 1.5 for the dense lines.
The qualitative interpretation of the magnetotelluric data presented in this paper is the analysis of selected transfer function invariants. In the magnetotelluric method, the impedance tensor can be transformed into seven independent invariants [40]. We used three basic invariants in particular [41]. The first one is related to the difference between the impedance tensor Z counter-diagonal elements. Hence, the contrast of the principal components of the impedance remains the same, regardless of the rotation of the whole system. The second invariant is based on the sum of the impedance tensor’s diagonal elements, the complementary components. It does not change regardless of the rotation angle of the system. The last of the fundamental invariants is the determinant of the impedance tensor, which also assumes a constant value when the system is rotated. Using the invariant, the so-called effective impedance can be obtained, which is calculated as follows:
The effective impedance can be transformed into the effective resistivity, and such values are presented in Figure 5. Analyzing the distribution for several frequencies, we see that the conductivity distribution structure in the studied area is three-dimensional, although, for greater depths, the trend is two-dimensional. On the southwest side, a high-resistivity complex corresponding to the Sudetes block is revealed.
Similar conclusions can be drawn from the analysis of the invariants of the Horizontal Magnetic Tensor. The HMT is a relationship between the horizontal magnetic fields at the observation point and the reference station ([18, 42]; Egbert & Kelbert 2012, [43]). A very interesting feature of the HMT is its ability to map the location of well-conducting rocks using the spatial distribution of certain rotational invariants. The figure (Figure 6.) displays the most informative HMT invariant (i.e., the largest singular values), which corresponds to the maximum values of the induced magnetic field. Larger values of this invariant correspond to a greater concentration of telluric currents, corresponding to the conducting structures.
However, the image is slightly distorted by the currents “flowing around” the nonconductive complexes. Therefore, the image is more complicated for smaller depths because we have complex resistivity distribution in the subsurface layers. Deeper down, the situation is stabilizing, and it can be said that the structure is 2D with local 3D anomalies. Here, too, the nonconductive block of the Sudetes is visible.
6. Results and Discussion
The results of 1D and 2D modeling will be compared with the full 3D model presented in [44]. A parallel variant of the latest ModEM 3D inversion code [42] was used, which employs MPI (Message Passing Interface) and which, besides impedances, includes tippers (Egbert & Kelbert, 2012). It uses a three-dimensional variant of the NLCG method [43]. A computing cluster with a standard MPI belonging to the Institute of Geophysics PAS was used for calculations. It was assumed that the covariance matrix values in horizontal directions were 0.7 but decreased gradientally to 0.5. In a vertical order, the value of the parameter was set at 0.7. As in 2D cases, the starting model of a homogeneous hemisphere with a resistivity of 100 Ωm was used. Such an approach made it possible to obtain a three-dimensional electrical resistivity distribution model, which corresponded to geological and geophysical assumptions in the most realistic way. As a result, a full 3D model of the conductivity distribution in the studied area was obtained (Figure 7).
To carry out the three-dimensional inversion, 49 of 51 MT stations were used. Just like during the two-dimensional inversion, the two most distorted stations were rejected. For each of them, 12 periods were selected, four of which fell within one logarithmic decade. The central area where all soundings are contained is 150 by 225 kilometers in horizontal directions with a mesh of 3 kilometers. The net’s total size with edges was 1350 by 1425 kilometers and 600 kilometers vertically. The total number of meshes was 68 out of 93 out of 47 in appropriate directions, with the first ten meshes in the vertical direction were in the air layer, which is hidden in the program code. Inverse calculations were completed after 75 iterations, with the final RMS matching error of 2.34. Lower impedance tensor errors of 7.5% and more minor errors for tippers of 0.015 are assumed.
Figure 7 represents the result of the three-dimensional inversion with main tectonic lines and the location of measuring points. A more detailed geological analysis of the results obtained from this inversion is described in another paper [44].
The dominant feature of the obtained resistivity distribution model is well-conductive structures in the northeastern part of the studied area and highly resistive in its southwestern part. The results confirm the location of the Dolsk and Odra faults postulated based on geological data, and they also show that these are very deep faults reaching the lower crust.
For both 1D and 2D inversion, an essential part of this step is the soundings at an appropriate angle. The -axis should hit the central axis of the structure under study. In this case, all soundings were rotated by an angle of 42 degrees. This angle was determined based on geological data and corresponds to the main fall of the two investigated structures. Moreover, the analysis of the polar diagrams and phase tensor also confirms this arrangement of geoelectric structures. It is dominant for most of the soundings, however, not for all. However, a constant rotation angle is assumed for the entire profile. This rotation allowed the TE and TM to be determined accordingly. For better comparison, the 1D inversion results are presented in section form.
The inversion was made on five parallel profiles of different lengths (Figure 1). The shortest profile marked as P0 is 123 km long and contains eight magnetotelluric soundings. It is also the most northern profile. The profile numbering increases to the southeast, while the estimate of the soundings to the southwest. The most extreme profile is marked as 0 because it was made additionally and is much shorter than the others. The P1 profile is 175 kilometers long and consists of nine soundings, another profile (P2) is slightly longer at 181 kilometers with the same number of measuring points. The most extended profile is the one marked as P3 and is about 187 kilometers long. It contains ten magnetotelluric soundings equally. The last profile (P4) is 180 kilometers long, with twelve soundings. Maps of all profiles have been made on the same scale and colour palette, ranging from 1 to 10000 Ωm. The same range was used to present the results of two-dimensional and three-dimensional inversion. The elevation was marked on each of the sections.
The resistivity distribution on pseudo-2D cross-sections for 1D smooth inversion (Figure 8) clearly shows two zones with significantly different resistivities. There are shallow conductive areas in the northeastern parts of almost all soundings, marked in red. It is not visible on the P0 profile, but it may be shorter and not cover the described zone. There is a significant resistivity zone in the middle part, clearly visible on the P1 and P3, and slightly less on P2 and P0. It is not visible on the last of the discussed profiles. This area is probably connected with the Wolsztyn-Leszno High (WLH), constructed mainly from poorly conductive rocks. The distinction between an area with anomalously high resistivity and a well-conducting block in the northeastern part of the studied area coincides with the assumed occurrence of the Dolsk fault, is also clearly visible. A horizontal layer of slightly higher resistivity is visible on practically all profiles in the well-conducting block described earlier at the northeastern end of the area. A considerable geoelectric differentiation of structures is evident on each profile but must continue on adjacent profiles in a similar location.
In the last two profiles, a high resistive background is visible at a depth of about 30 to 40 kilometers in the central (P3) and southwestern part of the profile (P2 and P4), which sinks into the northeast. Moreover, apart from the P3 profile, the high resistive background seems to be visible at the same depth in the northeastern part of the studied area.
There are two layers of high resistivity and one of markedly low resistivity. Among the higher resistivity structures, the substrate at about a 30-kilometer depth is the most pronounced, this time along virtually the entire length of all profiles. The low resistivity layer, located in the middle, is most clearly visible between about 10 and 20 kilometers, focusing on maximum conductivity in the final part of the profile, where it reaches the surface and significantly increases in thickness. The resistivity distribution coincides very well with the mapped geologic boundaries.
For quantitative two-dimensional interpretation, a smooth inversion algorithm—NLCG (Nonlinear Conjugate Gradient)—was used. A detailed description of the method, smoothing parameters, and the number of iterations is described in the previous chapter. The two-dimensional inversion was performed on the same five parallel profiles as the one-dimensional inversion (Figure 9). The lines are numbered from northwest to southeast, i.e., precisely the same way as before. The first two profiles, P0 and P1, clearly show two zones with significantly different resistivities. The first is the surface low-resistive structure with the minor resistivity focus in the profiles’ northwestern part. The second is a deeper structure with increased resistivity, with the highest saturation in the middle of the profile, at a depth of about 10 kilometers. It corresponds to the results obtained with Occam’s one-dimensional inversion. It is most probably connected with the Wolsztyn-Leszno block’s presence, built mainly of poorly conductive rocks.
The previous two figures (Figures 8 and 9)show the 1D and 2D inversion results obtained with WingLink and presented graphically as in the software. The following figures give the same results but are presented slightly differently and with a different colour scale. That is why it would be possible to compare these results with the results of 3D inversion. Figure 10 shows the presentation of 1D Occam inversion results in the depth range from 0 to 50 km on individual profiles.
Despite a different colour scale and a slightly different spatial arrangement, the division of the area between the well-conducting NE and the high-resistive SW is clearly visible. Moreover, an apparent boundary between the high-resistance SW block, a structure with intermediate resistance in the central part, corresponds to the arrangement of the Odra fault.
Figure 11 presents 2D NLCG inversion results in the depth range from 0 to 50 km on individual profiles, in the same scale and arrangement as the previous 1D result.
Similarly, in the case of two-dimensional inversion, low resistivity areas (high conductivity) are marked in red in the northeastern parts of almost all profiles. The distinctive closing of structures into the circles’ shape is associated with a considerable distance between the individual soundings. Complete mitigation would be related to a very high coefficient of horizontal smoothing, which could introduce some distortion and not reflect the area’s actual geological structure. A substantial geoelectric differentiation of structures is visible on the first two profiles, and larger systems are continued on the other lines in a similar location. Sharp delineation of structures is not well observed in this case or at least not as well as in the results of one-dimensional inversion.
When analyzing deeper sections, it seems surprising that there is no stable resistor in the background expected after analyzing 1D inversion results. It is slightly marked on the P1 and P3 profile, while on the P0 profile in the southwestern part, a structure of good conductivity is visible at a depth of about 60 km.
The best current method to achieve a complementary interpretation in the magnetotelluric method is three-dimensional modelling. The complex geological structure of the studied medium is dictated by the domination of three-dimensional structures, which is confirmed by high values of obliqueness and very polar diagrams, described in detail in previous chapters. We did not use the results of the lower inversions as starting models for the higher inversions, and our goal was to compare the results for each inversion separately.
In the initial phase of modelling, several preliminary tests were carried out to find the most reliable inverted three-dimensional distribution model of the electrical conductivity. A freely available ModEM software package by Egbert and Kelbert was used. A homogeneous half-space with a resistivity of 100 Ωm was used as a starting model. Different starting models were tested, which is very important to carry out correct three-dimensional inversion. A detailed description of 3D modeling can be found in [44].
The previous figure (Figure 12) shows the distribution of resistivity in the obtained 3D model, presented in the form of cross-sections through individual profiles in the depth range from 0 to 50 km. On the above distribution, very similar structures are visible to those obtained for one- and two-dimensional inversions. However, a much better fit for geological formations has been observed, and it is almost the same as the Odra fault zone. The P2 and P3 profiles also show a solid, weakly conductive structure, bordering on the northeast with a strongly contrasted, well-conducting area. This low-resistive area visible on the three central profiles coincides with the assumed location of the Dolsk fault.
Similarly, a well-conducting structure in the northeastern part of the study area is well visible. It is situated in the outer part of the model, and boundary effects may have influenced its shape. However, the presence of reduced resistivity on practically every section included in this paper allows us to presume that this is connected with the Szczecin-Miechow Synclinorium, limiting the Fore-Sudetic Monocline from the north. The horizontal arrangement of most of the structures, especially in the initial part of the profiles, is also noticeable in the sections. The angle of the collapse of the layers reflects their homoclinic arrangement.
Comparing the results of 2D and 3D inversion, in the most suitable form, i.e., profile cuts (Figure 10 for 2D and Figure 12 for 3D), a certain correlation can be observed. Although we have a clear resistivity differentiation between the southwestern and northeastern parts in both cases, it is more pronounced in the 3D results. This is probably due to the geoelectric three-dimensionality of the structures studied and that the rotation direction chosen was only roughly hit. Nevertheless, if we do not have the computational power or time to perform a 3D inversion, the results of the 2D inversion give a rough approximation of the studied structure.
For comparison with the results of the other inversions, the modelling with the Occam algorithm for the effective impedance (invariant) is presented as the result of the 1D inversion in both cases (Figures 8 and 10). This means that the results obtained in this module should be compared with those obtained for each polarization separately in order to draw the correct conclusions. This is a perfect solution when dealing with highly distorted data. The depth slices (Figure 13) show slight variation between the first five kilometers. The picture is slightly different for 13 kilometers, and significant changes begin for much deeper layers. For a comparison with the 3D inversion results, the 1D inversion results will be presented here as depth slices from the same depths as 3D (Figure 14). Due to the nature of the 2D inversion and the modelling along with the profile, it does not make sense to rearrange the results as depth cuts for comparison with the other inversions. In this case, it was decided to compile the results of the 2D inversion arranged in the same way and at the same scale as the results for the 3D inversion (Figure 12).
Figure 14 shows three-dimensional models of resistivity distribution in the form of depth maps for selected depths. It is suitable for shallower surface layers and deep structures, and only those structures should be analyzed, which are located within the area covered by the soundings. The model calculated outside this area is subject to a substantial error. From a depth of about two kilometers, low-resistive structures, reflecting surface layers of the overburden with good electrical conductivity, begin to dominate. The high resistivity zone in the southern part of the studied area coincides with the black marked Odra fault. Below three kilometers deep, the well-conducting region begins to blur, limited to small clusters, mainly in the northeastern part of the study area. However, the areas with increased and high resistivity were also visible on the cross-sections and the depth slices for 1D inversion (Figure 9). Domination of highly resistive structures in the central part of the studied area at a depth of about 5 kilometers may be related to the Wolsztyn-Leszno High, quoted many times before when describing the results of other inversions. An important observation is that the high resistive structures do not exceed the line determining the Dolsk fault location, which may prove its actual occurrence in this area.
Among the deeper structures, less than 20 kilometers, a division seems to be noticeable between a high resistive area in the southwestern part and a well-conducting northeastern part. This division was evident basically for the results of all inversions carried out. Thus, a clear division into two highly different geoelectric centers in the deep parts of the studied area is visible.
When comparing depth slices of 3D (Figure 14) and 1D (Figure 13) models, more similarities can be found for approximate depths than in the analogous situation when comparing 2D and 3D models. In both cases, the spatial continuity of the structures is preserved, although, as before, much better rendered in the 3D model. In the 3D inversion results, regional structures such as the Odra fault or the Dolsk fault are much better indicated. Based on only a few soundings, the WLH marks better in the 1D model.
The two-dimensional inversion results give the least reliable and readable result, and instead, the 3D inversion results are the best. The 1D inversion, on the other hand, gave satisfactory results, somewhat close to the final 3D inversion results. Such a conclusion is possible because 2D inversion strongly depends on the arrangement of the structures and how oriented the profile is, which is subjected to the inversion.
7. Conclusions
The electrical resistivity of the Earth is the parameter that provides valuable information that distinguishes geological structures with different petrophysical characteristics. Such geophysical identification of structural boundaries is important for the studied area’s comprehensive geological-tectonic analysis and may translate directly into the verification of proposed scenarios of geodynamic evolution.
Comparing the modeling results with algorithms of different dimensions shows that the essential features of the geological structure are revealed in each case. However, the approximations of the extent and depth of these structures may differ significantly. In the case of smaller structures, the differences are even greater. Depending on the dimensionality of the algorithm, they may or may not be present in the final model, and their location and extent also depends on the algorithm used.
It is noteworthy that in geological conditions where three-dimensional structures dominate, one-dimensional inversion may give better results than two-dimensional inversion. This is because, in the case of two-dimensional inversion, to carry out a correct interpretation, it is imperative to locate the profile across the span of the studied structure. However, if there are more tested structures arranged at different angles or the measured form does not have the axis of homogeneity, the obtained results of two-dimensional inversion will be ambiguous and burdened with a significant error. It was the case for the area under study. For the correct interpretation of the MT data, 2D or 3D inversions should be avoided as the only approach to MT data imaging. 1D inversion is much needed to ensure the quality of 3D interpretation in combination with 1D as well as 2D and correct defining of the geological concept. The comparison of different inversion methods shows that 1D and 3D are closely related. The inversion results are comparable, except that one-dimensional inversion reproduces the high resistive structures but smears them, while 3D inversion sharpens the image considerably. It is challenging to separate specific depths and thicknesses for the results of single-dimensional inversion.
Regardless of the visible differences between the results of modeling with 1D, 2D, and 3D algorithms, they all confirmed the location of the Dolsk and Odra faults. The elevation of Wolsztyn-Leszno, strongly connected with the Dolsk fault, is also well observed. However, the elongation of the East European Craton to the southeast of the TTZ and its continuation in the Palaeozoic Platform’s deep background remain unclear. The results presented here do not indicate how far the recent lowering of the Baltica under cover of younger sediments reaches. However, the results obtained correspond to the current knowledge on the subject, which allows us to suppose that it concludes with a probability that borders on the Dolsk fault and perhaps even on the Odra fault [17, 45].
Data Availability
The edi data used to support the findings of this study are currently under embargo while the research findings are commercialized. Requests for data, 12 months after publication of this article, will be considered by the corresponding author.
Conflicts of Interest
The authors declare no conflict of interest.
Authors’ Contributions
“SO,” “WJ,” and “WK” were in charge of the conceptualization; “SO,” “WJ,” and “KN” were in charge of the methodology; “KN” was in charge of the software; “WJ,” “WK,” “KN,” and “WK” were in charge of the validation; “SO” and “WJ” were in charge of the formal analysis; “SO” was in charge of the investigation; “SO” was in charge of the resources; “SO” was in charge of data curation; “SO,” “WJ,” and “KN” were in charge of the original draft preparation; “SO,” “WJ,” “KN,” and “WK” were in charge of the review and editing of the manuscript; “SO” and “KN” were in charge of the visualization; “WJ” and “WK” were in charge of supervision; “WJ” was in charge of project administration; and “WJ” was in charge of funding acquisition. All authors have read and agreed to the published version of the manuscript.
Acknowledgments
This work was partially supported by the National Science Centre Grant (NCN) number 6702017/25/B/ST10/01348.