1. Introduction
The tectonic structure of the Antarctic continent is the least studied compared to the other continents on the Earth. The reason for this lies not only in the perplexing geography and the extremities of its climate, but also because it is covered by ice sheets with a thickness of over 2 km in a vast region [
1]. Thus, applying geophysical methods to explore subglacial structures could be an important tool to better understand the tectonic structure in Antarctica, for example, in addressing issues such as subglacial hydrologic systems and tectonic deformations. The most favorable methods allow large-scale and efficient data collection and smooth operations on the ice sheets without costly logistical support requirements. The magnetotelluric (MT) method is one of these approaches. It explores the subsurface electrical structure by measuring the induced electromagnetic signals on the surface that derive from the induction of natural alternating electromagnetic fields in the Earth. MT can help in subsurface exploration to the depth range from tens of meters to hundreds of kilometers depending on different periods of signals.
In addition, MT data are sensitive to temperature, fluids, and molten material within rocks. These features are closely relevant to the dynamics of the ice stream and the rheology of the subglacial rocks. Therefore, the MT method has been widely applied to study the water distribution below the thick ice cap [
2,
3,
4], deformations of the central Transantarctic Mountains [
5], and the magma activities beneath the Deception Island [
6] in Antarctica. These electrical resistivity structures derived from MT data have been playing an important role in updating our knowledge on issues ranging from large-scale tectonic deformations to causes of individual ice streams [
7].
The southern Prydz Bay region in Eastern Antarctica probably includes essential clues for amalgamation or/and splitting of Gondwana [
8]. In addition, it is a landfast sea-ice concentrated area [
9]. This region is one of the excellent locations for studying the mechanisms for continent evolutionary processes, ice flow, and mass loss. Undoubtedly, the MT method is a potential technique to provide new evidence for answering these questions. However, so far, very few MT surveys have been conducted in this area. The first MT survey, including two broadband measurements in southern Prydz Bay, was conducted by Kong et al., in 1994 [
10] in which a one-dimensional electrical resistivity structure was revealed beneath two sites. About twenty years later, Peacock and Selway [
11] reported a 3D resistivity model by inversion of MT profile data, which parallels the southern Prydz Bay coastline from the Vestfold Hills to the west of the Rauer Group area. An extensive low resistivity body was imaged in the mid-lower crust below the Rauer Group, which coincides with the high positive magnetic intensity observed on the surface. Based on supplements from the Chinese Zhongshan Station, eight long-period MT data points along a profile were collected in the south part of Larsemann Hill [
12] from 15 November 2017, to 5 February 2018, in China’s 34th Antarctic expedition. They interpreted these data by 1D Occam inversion but proposed that the data had shown the distinct 3D effect of the crustal resistivity structure.
A new broadband MT data set was collected at 11 sites during China’s 36th Antarctic expedition in 2020 (
Figure 1). Guo et al., modeled these data with 2D inversion and found conductors in the crust [
13]. According to a detailed analysis of the impedance tensor and tipper data, we recorded a clear 3D structure below the MT profile, for which a 3D interpretation was required. In this paper, a new 3D resistivity model derived from the 3D inversion of these data is shown. In addition, before the 3D inversion, we reprocessed the time series. The magnetic inclination of the data was corrected by using calculation results from the latest IGRF model [
14]. In order to obtain robust and high-quality MT impedance tensors and tipper data, we ran the EMTF code package [
15] to convert the time series to MT transfer functions. Meanwhile, the most reasonable curves were selected by comparing different results from different coherence sorting parameters. Besides, some other important aspects of the MT data are addressed in this paper.
2. Data Acquisition, Time-Series Processing, and Qualitative Analysis of MT Transfer Functions
A broadband (0.003125 s~1000 s) MT data set was collected along a profile at 11 sites from 12 December 2019 to 15 February 2020 (UTC). All these data were collected by an MTU-5A instrument produced by Phoenix Geophysics. The variations of the electrical and magnetic fields were recorded by nonpolarizing Cu-CuSO4 electrodes with 100 m electrodes spacing and by broadband MT coils Phoenix MTC-50, respectively. The occupation time on each site ranges from 25 h to over 6 days (the specific time for each station is shown in
Figure S1). This profile is over 20 km long extending to the south from Larsemann Hill (
Figure 1). Sites at the south end of the profile had been shifted eastward because all these sites have to be installed along a path that is confirmed safe to drive above.
Several authors [
7,
11,
16] have reported that MT data collection on the Antarctic continent probably needs three important problems to be taken into consideration. They are (1) the large contact resistance between the electrodes and the ice/slow ground, (2) the electromagnetic noises produced by the moving charged snow and ice particles induced by wind, and (3) the nonplanar magnetic source field, which conflicts with the fundamental assumption in the MT method.
In this paper, we considered these three problems in the acquisition and processing of the time series. First, we filled saturated saltwater into the electrodes’ hole to lower down the contact resistance between electrodes and slow ground (The measured contact resistance at most of the sites is below 1 MΩ). Second, the EMTF code package [
15], including robust estimation of MT transfer functions were applied to reprocess the time series, and the most reliable and clean curves were selected among different results from different coherence sorting parameters. In this processing, the effect of nonplanar sources was suppressed to some extent. In addition, the magnetic declination at each site was corrected by using calculation results from the latest IGRF model [
14]. This is very important for the MT measurement at the high latitude and polar areas, such as Antarctica, where relatively big magnetic declinations take place. Third, we visually estimated the apparent resistivity and phase curves calculated from the MT impedance tensors and looked at how these curves responded to the above-mentioned problems. We found that the sites were more or less affected by the noises at a period ranging from 0.1 s to 2 s, probably owing to the moving charged snow and ice particles. In addition, most of the data show an abnormal decline in the apparent resistivity curves at periods shorter than 0.01 s. It is a typical characteristic related to the high contact resistance between the electrodes and ice/snow ground [
17]. Nevertheless, these effects were small in the majority of period points, and we deleted those outliers induced by noises and omitted data points at short periods distorted by the high contact resistance. With regard to the effects from the nonplanar magnetic source, only the data points with a period shorter than 1000 s were kept for further analysis in which the nonplanar source effect was not present or largely suppressed by the robust estimation in the time-series processing [
18]. Besides, data measured at the most north (
Figure 1) were entirely removed from the data set due to the fact that they were over five orders of apparent resistivity higher than other sites (
Figure S2), which indicates this data were seriously distorted, probably by the small inhomogeneous bodies on the surface or by the seawater to the north. MT data at 10 sites were appropriate to be applied for qualitative and inversion. All the data are shown in the
Figure S2 in the Supplementary file.
A qualitative analysis of MT data was undertaken before the MT data were inverted. The results from phase tensor analysis are shown in
Figure 2 [
19,
20]. According to the variations in the axis directions of the ellipses at different sites and periods, the results indicate that the regional electrical resistivity structure is largely 3D. A 3D inversion is required to interpret this broadband MT data set. Coordinate invariants of the phase tensor are shown in ellipses with different colors, which imply distinct variations in subsurface resistivity at vertical and lateral directions. The data at periods 0.06 s indicate that the upper crust of the study area to the north is characterized by higher resistivities than the south segment. At periods 40 s and 125 s, the presence of higher values in the phase tensor implies low resistivities in the middle to lower crust. According to a longer period at 500 s, the decreasing values in phase tensor as compared to the periods at 125 s indicate relatively higher resistivities at depth, where probably the lower crust is present. These features can also be seen in the pseudo-section of the phase tensor in
Figure 2b, which shows that the study region generally consists of a resistive layer in the upper and deep lower crust, and a relatively low resistivity layer is sandwiched between them. In addition, the results of phase tensor in
Figure 2 indicate lateral variations in the mid-crustal low resistivity layer occur along the N-S direction.
The induction vectors derived from the vertical magnetic field are another critical tool to qualitatively analyze the lateral variations of subsurface resistivity structure. It is significant that most of the vectors, in Parkinson convection [
21], tend to point to the east-southeast at all the periods shown in
Figure 2a. This observation is important and has implied the following three facts: (1) The geoelectrical strike in the deep crust is nearly North-South, which corresponds to one of the directions shown in the rose diagram in
Figure 1. The rose diagram is a statistical result of the principal axis angles of the impedance tensor at different sites and periods. This means that the geoelectrical strike is almost subparallel to the profile and implies that a 3D interpretation is required for these data, (2) The seawater to the north has few impacts on the MT data at periods longer than 0.06 s. (3) There are low resistivities buried in the deep crust at the eastern side of the MT profile, as the induction vectors in Parkinson convention point to the low resistivities.
3. 3D Inversion
Inversion is required to convert the frequency domain MT data into a 3D model of subsurface resistivity. The MT data were inverted using the ModEM inversion [
22,
23]. The inversion model had a horizontal grid size of 1.0 km in the core area, which was covered by the data. The data selected for inversion included the full impedance tensor and tippers in the period range 0.003125 s–1000 s.
Grids of the inversion model in the vertical direction were assigned a thickness of 20 m for the first layer and increased to the depth of over 800 km with a factor of 1.15. The topography was not included in this 3D model as the surface is flat, and the seawater has less effect on the MT data according to the qualitative analysis described in the above section. There are 80 layers in total in the 3D inversion model including 10 layers above the surface set to model the air layers in ModEM inversion. In the horizontal direction, all the grids are parallel to the north-south or east-west. The cell size within the core area was set as 1.0 km × 1.0 km. There were 13 padding cells in four horizontal directions, with a constant padding rate of 1.5. The edge of the model was extended to over 600 km from the center, which met the forward boundary conditions of 3D MT inversion. In total, there were 54 cells in the north-south direction and 35 cells in the east-west direction in the model.
In the 3D inversion of these MT data, different data components, regularization parameters, and error floors were applied to invert the subsurface resistivity structure. According to these inversions, we found a small number of data points with severe distortions, which could not be fitted in all the inversion models and were therefore deleted.
In this paper, a preferred 3D model was obtained by 3D inversion of full impedance and tippers from an initial uniform half-space model of 3000 Ωm. This initial model’s resistivity value is chosen by considering the average level of the apparent resistivities at the shortest periods in the data. The model responses fit the MT impedance tensors and the trend variation of tipper data quite well (
Figure 3 and
Figure S3); even the absolute amplitudes of the tipper data are difficult to be fitted in the 3D inversion. In this inversion, the data error floors assigned to the diagonal and off-diagonal components were 5% of √(|Z
xy·Z
yx|), and an absolute value of 0.05 was assigned to the error floors of the tipper data. In addition, the smooth parameters of the three directions (
x,
y,
z) in the regularization term were set as the same value 0.3. We forced the inversion to update the regularization parameter when the relative change in the overall root mean square misfit (r.m.s. misfit) between two consecutive iterations was below 0.015%. The new regularization parameter was calculated automatically by using the last value divided by 10. The 3D inversion would terminate upon the regularization parameter decreasing to 10
−8, or when the overall r.m.s. misfit was below 1.05. Finally, this inversion began with an r.m.s. misfit of 36.74 and was reduced to 3.62 after 133 iterations.
4. Results
The preferred 3D resistivity model is shown in six horizontal slices at different depths and two vertical cross sections in
Figure 4 and
Figure 5, respectively. In general, the new 3D resistivity model reveals that a resistive upper crust overlays on a relatively low resistivity middle crust. This feature is consistent with the above-mentioned qualitative analysis results. Specifically, the 3D model shows a significantly low resistivity body (C1) extended in the crust, which migrated eastward as the depth increased. Even though the location of the C1 at depth is out of the MT profile moving to the eastern side, the 3D inversion is capable of constraining the nearby off-profile structures [
24]. In addition, the location of the C1 is consistent with the eastward-pointed induction vectors and high values on phase tensors at periods longer than 1 s (
Figure 2). In addition, C2 to the south at the shallow part (<6 km) also coincides with the high phase tensor values (Φ
min > 45°) at period 0.06 s as shown in
Figure 2a. Furthermore, the south end of the vertical north-south cross section shows a clear variational pattern in resistivity that is high-low-high-low resistivities in the crust column from the top surface to the depth (
Figure 5). This pattern coincides with the variation in the phase tensors of three MT sites at the south end of the profile as shown in
Figure 2. Thereby, this preferred 3D model has recovered the main features in the MT data required by the data. Besides, by comparing our 3D resistivity model with the seismic velocity model reported in previous studies [
25], we found that C3 in the shallow part (~0.2 km to 1 km) is consistent with a low-velocity zone revealed by the shear-wave velocity model that was constrained by the Rayleigh wave dispersion curves. The information is demonstrated in detail in
Figures S4 and S5. We can, therefore, conclude that the features in this 3D electrical resistivity model are robust and reliable. It is noted that there are low resistivities near the north and south ends of the profile. These structures are weakly constrained by the data, and are not interpreted in this paper.
In the previous 2D interpreted model, extensive low resistivities were shown in the middle to low crust to the north of this profile [
13]. After carefully analyzing this MT data set in this paper, we found that these data were seriously impacted by 3D structures in the subsurface, and the geoelectrical strike subparalleled the direction of this profile. Thereby, artificial structures will very likely be placed in the 2D interpreted model. In comparison, the 3D interpreted model was rather more robust and contained few artificial structures.
5. Model Interpretations and Discussions
C1 is the main feature revealed by this new 3D resistivity model, which extends eastward in the entire crust from the top surface to the depth at low crust (>20 km) as shown in
Figure 5. This feature has been discussed in detail below in this paper.
The C2 has its upper interface at depth of ~2 km. It most likely connects with the C3 at a shallower depth. Hill reported similar features revealed by a couple of MT data set on Antarctica [
7]. These low resistivities were attempted to be explained as interconnected fluids near the base of the ice cap [
5,
11,
16] because their depth coincides with the ice thickness (~2 km) within most of Eastern Antarctica [
1]. In contrast, the thickness of ice sheets in our study area was reported to be lesser than 700 m [
26]. In addition, these low resistivities coincide with a low shear-wave velocity zone (
Figure S5), and Fu et al., interpreted this shallow anomaly as intrusion rocks associated with Pan-African thermal events according to their velocity model [
25]. However, the dry-aged granitic rocks are usually characterized by over 1 MΩm in the absence of interconnected partial melts or pore fluids [
27]. Since we couldn’t find any evidence of thermal activities in later Cenozoic in this region and a low heat flux at present day was reported by Martos et al. [
28], we inferred that C2 and C3 in the shallow part (<6 km) of the crust were more likely caused by the water generated at the base of the ice cap that migrated downward into the depth through fractures in the shallow crust. In addition, the fluid fraction of 0.004% in an intrusive granitic rock would produce 1000 Ωm for the bulk resistivity [
29] (Park and Mackie, 2000), which matches the resistivity range of 300 Ωm~1000 Ωm of the C2 and C3.
The location of the dominant feature C1 coincides with the high positive magnetic intensity belt named Amery Lineament [
30] shown in
Figure 6. The extended direction of the C1 shown in
Figure 4,
Figure 5 and
Figure 6a indicates that C1 most likely sources from the east of the study region in the lower crust. Peacock and Selway measured MT data further eastwards of our study area [
11]. A significant low resistivity body coinciding with high positive magnetic intensity was discovered in the middle to lower crust beneath the Rauer Group. It related the low resistivities to the high positive magnetizations. In addition, according to the magnetic intensity map (
Figure 6b), a distinct lineament featured by high magnetizations named Amery Lineament is present in the southern Prydz Bay region, which runs from the southern Vestfold Hills across the northern Amery Ice Shelf and further west to the southern MacRobertson Land [
30,
31]. The MT profile in this paper crossed this east-west lineament, and a low resistivity C1 was revealed to extend to and connect with the east of the MT profile. Therefore, by integrating the study results of [
11] to the further east area and the magnetic intensity observations, we propose that the C1 in our model revealed the location and geometry of this linear high magnetic anomalous zone at depth when it laterally crosses the study area in an east-west direction (
Figure 6). Besides, as the thickness of the crust in this area is about 30 km [
32], we could infer that the Amery Lineament is probably at least a crustal-scale structure, and whether or not it has extended to the upper mantle is still unknown due to the low resolution of the data at depth over 30 km.
In the previous studies conducted on this area, the authors interpreted the low resistivities in their models as an implication of thermal activities supplied by heat sources from the mantle [
12]. However, a very low heat flux and a deep Curie depth in the study area were measured and calculated at the present day [
28], which does not support the interpretation of thermal activities in the crust. In addition, the hot intrusive igneous rocks resulting in fractures saturated by fluids are an alternative explanation for the low resistivities [
13]. In this area, however, most of the granitic rocks that outcropped on the surface could be dated at over 500 Ma [
33], and no evidence for extensive fluids having migrated into the crust below the depth of over 5 km was found. If any water was present in this area, the content is perhaps little and in the shallow part of the curst. To sum up, so far, it is more reasonable to associate the C1 with the high magnetic intensity observed on the surface, which could be because of sulfides, magnetite, or graphite in the granitic dikes in the deep crust [
11].
Finally, it is worth mentioning that our new 3D resistivity model revealed that the C1 outcropped on the shallow land surface was probably covered by an ice cap hundreds of meters thick. Shallow-drilled well data integrated with this 3D resistivity model will help in providing new insight on the compositions of the Amery Lineament in the future.