Validity of the dispersion relations in magnetotellurics. Part II: synthetic and field data

The impedance tensor and tipper are shown to be non-causal in some classes of 2-D and 3-D magnetotelluric (MT) models, including those widely encountered off shore. At the same time, examination of a large database of onshore MT records yields only a handful of non-causal tensors, being non-causal due to the influence of nearby artificial conductors rather than some exotic geological conditions. This result implies that the actual chance to encounter a non-causal impedance tensor or tipper in terrestrial MT exploration is vanishingly small, thus securing the universal applicability of the dispersion relation technique for quality assessment of the mentioned transfer functions.


Introduction
In the past century, the dispersion relations (DR) were widely used for processing, quality assessment and checking the consistency of broadband magnetotelluric (MT) data for qualitative interpretation (Kunetz 1972;Boehl et al. 1977;Fischer and Schnegg 1980;Sutarno and Vozoff 1991). However, in recent times, the technique has slowly started losing favor with the MT practitioners. The main reason is that the DR are universally valid only in 1-D (Weidelt 1972) and H-polarized 2-D (Weidelt and Kaikkonen 1994) isotropic media, while in more complex MT models with high resistivity contrast, they are sometimes found violated (Berdichevsky and Dmitriev 2008;Alekseev et al. 2009;Ichihara and Mogi 2009).
The arising problem is straightforward-is it legit to employ the DR violation criterion as a plain indicator of noisy/biased/poor processed data in a general 3-D case, if the complex geology alone can lead to DR violations? The main goal of our research is to substantiate positive answer to the above question by showing that any intrinsic (i.e., related to complex geology rather than noisy data) DR violation can be steadily identified. The other goal is to estimate the occurrence frequencies of the intrinsic DR violations in field MT data.
We have spent the first part of the paper , hereinafter referred to as Part I) considering the issue of the DR validity in MT transfer functions from the most general perspectives of the black box approach. For ease of further reading, the main conclusions are briefly summarized below (the angular frequency is hereinafter denoted by ω , while X stands for an arbitrary MT transfer function in the frequency domain-e.g., the impedance, tipper, etc.): • Any tensor X at a given observation site is either causal or non-causal, and this property is invariant to the coordinate system rotation. • Any component X ij of a causal tensor X is either minimum phase (MP) or non-MP, and this property is sensitive to the rotation procedure. • The DR of the first kind (DR-I) Open Access *Correspondence: shimizu@eri.u-tokyo.ac.jp 2 Earthquake Research Institute, The University of Tokyo, 1-1-1, Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan Full list of author information is available at the end of the article is valid in all components of a causal tensor X . • The DR of the second kind (DR-II) is valid only in MP components of a causal tensor X . • The DR-II violation in any non-MP component of a causal tensor X is strictly positive (for the harmonic time factor e iωt used in the paper) and consists of one or more 180° phase flips, known in the MT literature as the "phase rolling out of quadrant" (PROQ) phenomenon. • The DR-II violations in a non-causal tensor X consist of at least one negative 180° phase flip, but may also include some positive 180° phase flips on top as well.
With that, the exact number and frequency behavior of the positive flips are sensitive to the rotation procedure and vary from component to component, while those of the negative flips do not.
From practical standpoint, the obtained results imply that any component X ij of any MT tensor X in a given coordinate system can be relegated to one of the three general classes of the DR validity (Table 1). According to this classification, the DR application for quality assessment of causal MT tensors acquired in arbitrarily complex geological conditions is straightforward: if there (1a) ImX ij = π 2 · dReX ij d ln ω * B(ln ω) = DR ReX ij , (1b) B(ln ω) = 2 π 2 ln coth ln ω 2 , (2a) arg X ij = π 2 · d ln X ij d ln ω * B(ln ω) = DR ln X ij , (2b) B(ln ω) = 2 π 2 ln coth ln ω 2 is no PROQ phenomenon, then any observed DR-II or DR-I violation comes from poor data quality; if a positive phase flip is revealed, then one should use the DR-I only (or, if deemed more suitable, modified DR-II-see Eq. 7 in Part I). Non-causal MT tensors, if they do at all exist, can be steadily identified and equivalently represented via one or several causal tensors from the same response space.

Modeling
According to Table 1, an intrinsic DR-II violation in any MT transfer function is always associated with some "unexpected" positive or negative 180° phase flips (i.e., sign reversals), which naturally suggests looking for such irregularities in the synthetic models where some of the MT field components reverse their initial direction with frequency.

1-D model with anisotropy
In isotropic 1-D media the off-diagonal impedance components Z xy and Z yx are causal and minimum phase (Weidelt 1972; see also Berdichevsky and Pokhotelov 1997;Weidelt 2003 for the case of polarizable layers), while the diagonal components Z xx and Z yy are identically equal to zero. However, a completely different situation arises if some of the layers in a 1-D model are considered to be anisotropic: in this case, Z xx and Z yy = −Z xx become generally non-zero and may show numerous sign reversals (Marti 2014). The simplest model we found to reveal a non-MP magnetotelluric response is given in Fig. 1. It consists of two anisotropic layers with the same values of principal resistivities ρ x ′ = 1 Ωm , ρ y ′ = 100 Ωm and ρ z ′ = 100 Ωm , but with different strike angles α strike (azimuth between the x-axis and the corresponding principal direction x ′ of the resistivity tensor): − 45° for the upper layer, and + 45° for the lower one. The forward problem was solved in the MT1DAnisoFwd software (author D. Alekseev) based on the numerical algorithm developed by Kovacikova and Pek (2002). The Z yy response calculated for the upper layer's thickness h = 0.1 km, 1 km and 10 km is shown in Fig. 2 (note using Basokur's frequency normalization Ẑ n =Ẑ/ √ iωµ 0 , discussed in Part I, for the lower curves). The depicted MT data clearly shows that Z yy represents a causal but non-MP function (class 2), with valid DR-I and distinctly violated DR-II. This violation is associated with the reversal of the corresponding component of the induced electric current, which takes place at the periods when the relative influence of the basement layer on the measured response outweigh that of the upper layer.
Since in a two-layered medium the impedance phase values represent functions of h √ ω (Kovacikova and Pek 2002), a tenfold increase of h naturally "shifts" the whole Z yy phase curve along with all its irregularities by 2 decades of period rightwards. As a result, the DR-II violation becomes practically observable only if the upper layer's thickness corresponds to the measured frequency range. For h > 1 km and T < 1 s the influence of the basement layer on Ẑ is negligible and the DR-II in Z yy is effectively valid. Similarly, for h < 0.1 km and T > 10 s, the subsurface layer ceases to affect the frequency behavior of Ẑ , Fig. 1 Anisotropic 1-D model. Principal resistivities for both layers are ρ x ′ = 1 Ωm , ρ y ′ = 100 Ωm and ρ z ′ = 100 Ωm Fig. 2 The DR validity in the Z yy component (class 2) for various values of the upper layer's thickness h and the DR-II in Z yy turns out to be valid as well, albeit after the corresponding correction for the initial sign of the measured response. This means that MP or non-MP behavior is essentially a local (in frequency) rather than global attribute of a transfer function, and the same component of Ẑ may be relegated to different classes of the DR validity, depending on the particular period range. For example, in the considered 1-D model with h = 1 km (Fig. 2, middle column) Z yy is effectively MP (class 1) within the AMT period range (from 10 −4 s to 1 s), but turns out to be non-MP (class 2) when considered at MT periods (from 3 · 10 −3 s to 10 3 s).
The DR-II violation in Z yy defined as and calculated for various thicknesses h of the upper layer is shown in Fig. 3 (dotted lines). All curves have the same smooth shape, which can be easily identified (cf. Fig. 9a in Part 1) as that of the elementary phase lag θ 1 , given by This implies that, when considered on the complex Ω-plane, Z yy (Ω) has a single zero right on the negative half-axis of imaginary frequencies, and its location Ω 1 = −iz 1 may be steadily defined by fitting the obtained data with the corresponding theoretic curves (Fig. 3, solid lines). The resulting estimates of z 1 (namely 2.8 · 10 1 rad/s for h = 0.1 km, 3.1 · 10 −1 rad/s for h = 1 km and 2.9 · 10 −3 rad/s for h = 10 km) characterize the central frequencies of the DR-II violation in Z yy for the given h values and expectedly differ from each other by about two orders of magnitude.

2-D model with topography
As follows from numerous published results (see Marti 2014 for a review), incorporation of strongly anisotropic objects into a 2-D model may lead to PROQ not only in diagonal, but also in off-diagonal impedance components (3) ϕ yy = arg Z yy − DR ln Z yy , as well. However, all such examples apparently belong to the same class 2 of the DR validity as the Z yy component in the above 1-D model, with the DR-II violated by the amount of an elementary phase lag θ n and the DR-I being universally valid. Emergence of a non-causal (class 3) impedance component requires strong perturbations of horizontal magnetic field (see Part I), which could be associated with the phenomenon of excessive concentration of the induced currents in highly conductive elongated structures, usually referred to as "current channeling" (Simpson and Bahr 2005;Berdichevsky and Dmitriev 2008).
Since the rigorous proof of the DR validity in MT impedance on the surface of isotropic 2D Earth exists only for H-polarized models (Weidelt and Kaikkonen 1994), the efforts of many researchers within the last decades were directed to numerical investigations of the MT curves in the vicinity of highly conductive E-polarized 2-D bodies (Berdichevsky and Dmitriev 2008). To the best of our knowledge, none of the proposed models have revealed any observable DR violation. However, in some synthetic examples with huge resistivity contrast (Parker 2010;Selway et al. 2012) the impedance phase curve was found out to be slightly leaving its quadrant downwards and then immediately returning back to "normal" values from the same side (i.e., without any phase flips). Though such behavior does not yet imply the DR violation (see Zorin and Alekseev 2018), it clearly indicates the existence of strong magnetic field perturbations, which cannot be encountered on the surface of H-polarized media. A possible reason why these perturbations do not lead to the breakdown of causality is that in flat models considered by the above authors the induced currents are flowing exclusively below the observation sites. As a result, the horizontal component of the secondary magnetic field produced by a 2-D current channel of any intensity has the same sign as the primary magnetic field, and hence cannot force the observed response to reverse its polarity. In this situation we may yet expect the total Fig. 3 The DR-II violation in Z yy (dotted line), the best-fit elementary phase lag θ 1 (solid line) and the corresponding lower half-plane zero (open circle) of Z yy (Ω) for various values of h magnetic field to be reversed at an observation site located below the conductive body, e.g., due to steep 2-D topography, as that shown in Fig. 4.
To check this assumption, we assembled the given model of an "elevated 2-D conductor" on a 220 × 80 finite-element grid with minimum mesh size of 2.5 × 5 m, and applied to it the NWMT2D code (author A. Kaminsky). The obtained MT response reveals dramatic violations of DR-I and DR-II in the considered transfer functions (Fig. 5), thus showing that both Ẑ and Ŵ observed near the cliff are non-causal.
As follows from Table 1, one of the possible ways to assess quality of non-causal (class 3) transfer functions using the DR technique is to consider the corresponding inverse tensors: Ŷ instead of Ẑ , Ŷ n =Ŷ √ iωµ 0 instead of Ẑ n =Ẑ/ √ iωµ 0 , and so on. Nonsquare tipper matrix Ŵ for this purpose might be replaced by the pseudoinverse matrix Ŵ + defined by the following equation: The spectral components of Y yx and W + yz are plotted in Fig. 6. Though these functions represent exactly the same data as contained in Z xy = 1/Y yx and W zy = 1/W + yz , respectively, they both belong to the (lower) class 2 of the DR validity and hence are both causal.
The obtained results also reveal an interesting aspect of the DR violations in non-causal 2-D MT tensors, which at first glance may even seem to be self-contradictory. Indeed, on the one side, we see that the impedance Ẑ observed near the cliff is non-causal (Fig. 5); hence, all its components formally belong to the class 3 of the DR validity and both dispersion relations in them are meant to be violated. But at the same time, as shown by Weidelt and Kaikkonen (1994), in any isotropic 2-D model with arbitrary topography, the impedance tensor component corresponding to H-polarization of the medium must be minimum phase.
To solve this apparent contradiction we note that, for the model in hand, Eq. (13 in Part I) is simplified as follows: where the H-polarization magnetic response function h xx is MP, while the E-polarization magnetic response function h yy is not and hence have at least one lower half-plane zero. As a consequence, each of these zeros is translated from h yy right into detĥ = h xx h yy , thus becoming a non-causal singularity of the whole tensor Ẑ . However, h yy also appears as a multiplier within the Z yx component (Eq. 6), so that all its poles and zeros are h xx e xy 0 , Fig. 6 The DR validity in Y yx and W + yz components observed near the 2-D cliff mutually canceled out, and the transverse impedance Z yx turns out to be MP in the given coordinate system.
To have a closer look at this phenomenon we shall make use of the rotation procedure. Figure 7 shows the apparent resistivity ρ xy , impedance phase arg Z xy , DR-II violation curve φ xy = arg Z xy − DR ln Z xy and corresponding best-fit locations of the lower half-plane poles and zeros of the Z xy component rotated 0°, 15°, 30°, 45°, 60°, 75° and 90° clockwise, respectively. As expected, the locus of two thuswise discovered poles (viz., Ω 1 ≈ −8.6i rad/s and Ω 2 ≈ −91i rad/s) does not show any dependence on the rotation angle α , and the shape of the associated negative part of the total DR-II violation φ xy remains the same for all values of α . On the other hand, gradual rotation of the coordinate system forces the zeros of Z xy (Ω) to move along the complex frequency plane, and, as could be anticipated from the emerging negative cusp of lg ρ xy , at α higher than ~ 30° a pair of such zeros approaches the real-frequency axis from above. At α ≈ 39° these zeros cross the ω-axis at ω 1,2 ≈ ± 20 rad/s (± 3 Hz), turning for a moment the cusp of lg ρ xy into an infinite spike and introducing an additional positive 360° phase lag into the total DR-II violation curve. At α ≈ 56° the zeros meet at the imaginary frequency axis, split up, and start moving from each other in the opposite directions. Finally, at α = 90°, i.e., in the proper H-polarization mode, both zeros coincide with the corresponding poles, the positive and negative phase flips fully compensate each other, and the dispersion relations in Z xy become valid (slight mismatch of the compared curves at the very ends of the available frequency range is attributed to the inherent properties of the DR calculation-see Part I for details).
The above example shows that a component of a noncausal MT tensor may indeed be minimum phase-for that to happen, its lower half-plane poles have to be exactly overlapped by the corresponding zeros. However, such MP property is not the same as that in a causal tensor, since an infinitely small rotation of the coordinate system in any direction would lead to its failure. As a result, it seems reasonable to leave the general classification given in Table 1 intact and consider the MP components of a non-causal tensor as those belonging to the third class of the DR validity with the dispersion relations being violated by the amount of zero (i.e., θ n ≡ θ m ).

2-D model with bathymetry
While placing an MT station below a 2-D conductor requires rather exotic geological conditions on the land surface, this is evidently a common practice for the seafloor observations, which implies that non-causal impedance tensors and tippers be encountered in marine magnetotellurics on a regular basis. Indeed, the unique behavior of the offshore MT data in the coastal regions, which is "seldom, if ever, encountered on land" (Constable et al. 2009), was revealed decades ago and successfully simulated with the help of simple E-polarized 2-D models (e.g., Vanyan andPalshin 1990, White andHeinson 1994). Characterized by numerous negative phase flips, this behavior is reported to be a plain indicator of the DR-II violation in the offshore impedances (Alekseev et al. 2009, Kapinos andBrasse 2011;Kaufman et al. 2014). Surprisingly, none of the authors tried to check the validity of the DR-I in such data since this would have marked a decisive end to the hot discussion of the unsupported statement of Yee and Paulson (1988) about the general causality of Ẑ in real MT data (see Berdichevsky and Dmitriev 2008) by providing a simple and geologically reasonable counterexample. In this section, we fill this gap by pointing out the non-causal behavior of Ẑ in a simple E-polarized model of the coast effect.
The model consists of a conductive (0.3 Ωm) semiinfinite 2-D "ocean" located on top of a more resistive (100 Ωm) "continent", as shown in Fig. 8 (bottom). The origin of the horizontal axis y is chosen at the place where the ocean depth reaches half (0.5 km) of its maximum (1 km) value, which happens to be a convenient reference point for estimating the spatial limits of the coast effect influence on MT data (Worzewski et al. 2012). The forward MT responses were calculated using the NWMT2D code at 101 observation sites equally spaced along the 50-km survey line, with an increased period density of 15 per decade, necessary for accurate representation of the anomalous curves.
The map given in the middle part of Fig. 8 shows the absolute value of the DR-II violation in Z xy defined as This parameter is restricted to lie within the interval [0, π ] and thus is more convenient for examination of the MT responses with particularly complex frequency behavior, than the simple sign-variable parameter φ xy = arg Z xy − DR ln Z xy used above (cf. Figs. 3, 7). As seen from the map, the DR-II is prominently violated at all underwater stations from the very coastline and up to those located at about y ≈ 30 km. At the same time, the frequency behavior of this violation varies greatly with distance, involving almost all periods within the considered range of 1 s-10,000 s near the foot of the continental slope, but becoming sharply localized around the periods of 100 s-200 s at the remote offshore sites. To further investigate this phenomenon we have examined each survey site separately (7) |φ| xy = arccos cos arg Z xy − DR ln Z xy . Zorin et al. Earth, Planets and Space (2020)  The absolute values of the obtained complex frequencies at which Z xy (Ω) has a pole (or a pair of poles) are plotted with correspondingly labeled hatched lines over the DR-II violation map in Fig. 8 (middle). These curves are, in fact, describing how the non-causal poles of Z xy (Ω) change their locations in the complex frequency plane when we travel along the survey line (at any particular station their locations are securely fixed, like those of the zeros shown in Fig. 7). Now, going towards the coastline from the open water, we can see that at y ≈ 30 km, a conjugate pair of poles crosses the real frequency axis at ω ≈ ± 0.035 rad/s ( T ≈ ± 180 s) from above, resulting in the breakdown of Ẑ causality. At y ≈ 19 km, the poles reach the imaginary axis, split up, and start moving from each other, making the period range with the out-of-quadrant impedance phase wider. Around y ≈ 10 km, one of these poles finds yet another pole to combine with and cross the real frequency axis again at ω ≈ ± 5.2 rad/s ( T ≈ ± 1.2 s), while the other one stays alone up to the very coastline, where it finally gets out of the available period range towards higher frequencies.
Non-causal behavior of Z xy component at the MT stations marked with red color in Fig. 8(bottom) is also confirmed by Fig. 9, which shows the map of normalized DR-I violation parameter for Z xy defined as As follows from the figure, at periods corresponding to the non-causal poles of Z n xy (Ω) , the DR-I fails to correctly "predict" the actual sign of ImZ n xy , so that the inaccuracy of its application comes up to ~ 200% of the total Z n xy amplitude, which is in good agreement with the general expectations of the DR-I violation in this type of functions (see Appendix 2 in Part I). Figure 10 shows the ρ xy values along the survey line. As should be expected, the regions where the complex plane poles of Z xy (Ω) approach the ω-axis are characterized by emerging positive cusps on the resistivity curves, which turn into infinite spikes at the points where the poles are purely real, i.e., where the corresponding magnetic field component totally vanish. According to the papers (Key and Constable 2011;Worzewski et al. 2012), in simplistic 2-D models of the coast effect, one such point is always present at the distance r ≈ 0.09ρ l (h/ρ w ) from the middle of the continental slope and the period τ ≈ 0.19ρ l (h/ρ w ) 2 , where ρ w and ρ l stand for the water and land resistivities, respectively, and h is the ocean depth in kilometers. Substituting ρ w = 0.3 Ωm, ρ l = 100 Ωm, and h = 1 km into the above equations, we get the values very close to those describing the apparent resistivity peak observed at y ≈ 30 km. Though this particular singularity of ρ xy is not unique for the present model (we see another one at about y ≈ 10 km, T ≈ 1 s, and probably, there are even more of a kind at higher frequencies), it appears to be of most practical importance as the one drawing a clear boundary between the MT sites with causal and non-causal behavior of the impedance tensor at long periods.
In conclusion, it should be noted that the transverse (H-polarization) impedance value Z yx obtained for the given periods and model parameters demonstrates MP behavior at all MT stations along the survey line. This means that the same mechanism of "compensating zeros" encountered in the above model of an elevated 2-D conductor (Fig. 7) applies to the H-polarized seafloor models as well. Furthermore, the observed MP behavior of the horizontal electric field components ensures causality of the admittance tensor Ŷ = 1/Ẑ along the whole survey line, thus solving the issue of the DR application for quality assessment of the MT data heavily distorted by the 2-D coast affect.

3-D model with a curved conductor
For completeness we shall consider one of the curved conductive 3-D structures known for producing MT responses with PROQ (e.g., Weckmann et al. 2003;Thiel et al. 2009;Ichihara and Mogi 2009;Kaufman et al. 2014). The most simple and, hence, geologically reasonable among them is the L-form model of Ichihara and Mogi (2009). As shown by these authors, anomalous behavior of Ẑ over the L-shaped conductor is associated with local "reversed channeling" (sign reversal) of the electric current observed in a minor region of the model, while the horizontal magnetic field "varies slightly in all area". This result suggests that in the L-form model one can encounter the impedance components only of the first two classes of the DR validity, while the emergence of a non-causal impedance tensor on the flat Earth's surface requires more exotic shape of the conductor, such as the "S" shape proposed by Aleksanova and Blinova (Kaufman et al. 2014). Indeed, in the S-form model the reversed current is observed in a very large region and happens to be so intense that we should expect the associated secondary magnetic field to dominate the total magnetic field within this region, which may even cause its sign reversals along with the breakdown of Ẑ causality.
To observe the MT responses produced by both L-shaped and S-shaped conductors at once we have constructed a combination of these two models shown in Fig. 11. The resulting "SL-shaped" conductor (0.5 Ωm) is located 50 m below the flat Earth's surface, has constant thickness of 400 m and extends infinitely along both directions of the y-axis. The background medium is a resistive (500 Ωm) homogeneous half-space.
The forward problem was solved for 57 periods on a rectangular grid of 627,200 (112 × 112 × 50) meshes using the ModEM code of Egbert and Kelbert (2012). To reduce excessive time consumption of the corresponding calculations, they were performed on the supercomputing complex Lomonosov at Moscow State University (Sadovnichy et al. 2013). Figure 12a, b shows the maps of maximum absolute values of the DR-II violation (Eq. 7) in the off-diagonal impedance components Z xy and Z yx within the period range from 10 −3 s to 10 3 s. The same maps obtained for the coordinate system rotated 45° clockwise are shown in Fig. 12c, d. The validity of the DR-II in Ẑ components depends on the employed coordinate system, and this dependence also considerably varies from site to site. Particularly, we have discovered that for each point over the L-shaped part of the conductor (e.g., in the site A), there exist some rotation angles at which both Z xy and Z yx are minimum phase. In more heavily distorted area over the central part of the conductor (e.g., in the sites B and C) at least one of these components happens to be non-MP, irrespective of the rotation angle. Trying to delineate this area, we found out that a reasonable approximation could be obtained by highlighting the region of the DR-II violation in the Berdichevsky impedance Z brd (Fig. 12e), which is basically an arithmetic mean of the principal impedance components with an intrinsic property of being a rotational invariant (Berdichevsky and Dmitriev 2008): The DR-II violation map for another important rotational invariant-the magnetic tensor determinant detM -is given in Fig. 12f. As pointed out in Part I, both Ẑ and Ŵ are causal at all MT sites with MP behavior of detM obtained with the base site placed in an area with regular field behavior (i.e., over the 1-D background or, simply, far enough from strong magnetic anomalies). Consequently, the given DR-II violation map for detM implies that all components of Ẑ and Ŵ observed at the sites A and B should belong to either of the first two classes of the DR validity, while at the site C, we should encounter only those of the class 3, which is fully confirmed by Fig. 13.
For being able to apply the DR technique for quality assessment of non-causal MT data, it is necessary to use some other transfer functions, which contain the same information but are causal. In accordance with Table 1, for this purpose one may try to consider the corresponding inverse or inter-site transfer functions from the given (9) Z brd = Z xy − Z yx 2 . response space. The first approach, proved useful in the above 2-D examples, is not applicable to the present model, since all of the MT field components observed over the middle part of the SL-shaped conductor are greatly distorted. However, if the MT measurements are synchronously taken at some remote station with regular field behavior, it is always possible to express local MT transfer functions via the corresponding causal inter-site transfer functions. For instance, using the data from the base site shown in Fig. 12, the pair of tensors Ẑ and Ŵ at sites A, B and C can be equivalently represented by causal tensors M , Q =ẐM and Ŝ =ŴM (see Eqs. 12-17 in Part I), whose components belong either to class 1 or class 2 of the DR validity (Fig. 14).

Field data
As follows from the performed modeling, the field MT observations are potentially able to yield non-causal impedances and tippers. To assess the practical risks of encountering such data, the following search algorithm was applied to the available field data sets: 1. First, the impedance tensors with arg Z xy and/or arg Z yx leaving their regular quadrants are marked out. To make this procedure resistant to the possible presence of the outliers and other data inaccuracy, we use the criterion of the phase values being continuously out of quadrant for at least half decade of periods. 2. Then among the off-diagonal impedance components with out-of-quadrant phase values are chosen those with distinctly violated DR-II. As a robust indicator of such violation we use the condition p.v.φ ij > π/2 for at least half decade of periods, where p.v.φ ij stands for the principal value (i.e., one belonging to the interval from −π to π ) of the DR-II violation function φ ij = arg Z ij − DR ln Z ij . Thus- Fig. 13 The DR validity in Z yx and W zx components observed in MT sites A, B and C over the SL-form model wise discovered impedance tensors are labeled as "anomalous". 3. At the third step, all the determined anomalous tensors are further checked for validity of the DR-II in their rotational invariant Z brd (Eq. 9). The impedance tensors which satisfy the condition |p.v.φ brd | > π/2 on at least half decade of periods are thereby labeled as "strongly anomalous". 4. Finally, the normalized DR-I violation parameter Im n Z n brd (Eq. 8) is calculated for each strongly anomalous data example. All the revealed tensors with Im n Z n brd > 100% for at least half decade of periods are acknowledged to be non-causal.
The proposed algorithm is constructed in such a way that, being applied to the impedance tensors observed over the SL-shaped model (Fig. 12), it labels the sites A, B, C as anomalous, sites B, C-as strongly anomalous, and site C-as non-causal, while paying no attention to the (regular) Base site.

Land MT exploration
The overall MT database associated with numerous onshore projects performed by Nord-West Ltd. in the past 2 decades consists of more than 100,000 standalone AMT, MT and BBMT (AMT + MT) records acquired mostly in Asia, Europe and Latin America. However, for the purpose in hand we have removed from consideration relatively narrow-banded AMT records (so that only MT and BBMT ones are used in the below statistics) and repeated measurements acquired at the same location (so that for each unique survey site only the most recent record is used). At the same time, all the repeated measurements characterized by a considerable (> 50 m) shift of the survey location were left intact. Though such doing naturally introduces a number of low-quality records which then have to be identified and removed from the statistics in order not to be falsely labeled as anomalous, it ensures that none of truly anomalous MT responses be missed.
The result of applying the above algorithm to the whole selected database (60,226 impedance tensors, as of April 2019) and to some particular objects from it individually is given in Table 2. As to be expected, the overall percentage of the encountered anomalous impedances varies greatly from project to project and is also highly non-uniform in space (Fig. 15). In thick sedimentary basins and other "calm" geological conditions, MT surveying rarely yields more than 1% of anomalous data examples. However, as revealed by regional MT exploration projects considered here (e.g., the geotraverses 1-SB, 3-DV) and elsewhere (e.g., Weckmann et al. 2003;Thiel et al. 2009;Pina-Varas and Dentith 2018), this proportion may easily bounce up to more than 10-20% in orogenic belts, graphitized zones of ancient crystalline rock formations and other regions characterized by abrupt lateral variation of resistivity.
As follows from Table 2, pretty much every single impedance tensor (viz., 60,220 out of 60,226) from the database is labeled as causal, which implies that virtually all revealed MT data with PROQ belong to the second class of the DR validity. An example is given in Fig. 16. This record was acquired in August, 2017 at the site 28-002 of the Chara-Tokko MT project in Southeastern Russia (Fig. 15). Since only electric field was measured directly at the survey site, and the magnetic field taken from an adjacent location, the resulting transfer function formally represents an inter-site impedance Q (Eq. 16 in Part I). With phase rolling out upwards and DR-II violated by the presence of a ~ 180° lag at longer periods, the Q xy component is easily recognized as that belonging to the class 2 of the DR validity with a single non-MP zero ( n = 1 ). In Part I we propose three ways of quality assessment of such data using the DR technique: • apply the DR-I instead of DR-II; • modify the DR-II by the elementary phase lag θ 1 (Eq. 4); • rotate tensor Q , and if its components become MP, apply the DR-II to them.
The last approach happens to be of no use for the given case-the MT response observed at the site 28-002 belongs to the list of "strongly anomalous" impedances (Table 2), and either Q xy or Q yx is always non-MP, irrespective of the applied rotation angle. At the same time, application of the modified DR-II with an automatically fitted variable |Ω 1 | (Fig. 16, top) and straightforward DR-I application (Fig. 16, bottom) both work well. Now, let us consider the six non-causal impedances revealed in our database by the applied search algorithm. The first one (Fig. 17), characterized by much better quality than the other five, was acquired in October 2009 at the site 932G of the 3-DV geotraverse (Fig. 15), very close (< 100 m) to a local railway under construction. Since the railroad had not been put into service by the time of the field work commencement, it was not initially supposed to become an obstacle for MT exploration. However, the first obtained results proved just the opposite-apparently, in the context of highly resistive rock formations of that region, the secondary magnetic field associated with telluric currents flowing in the rails turned out to be intensive enough to overwhelm the primary magnetic field in near surroundings. As long as the initial location of the survey site happened to be on a hillside well below the railway embankment, the magnetic field distortion also resulted in the breakdown of causality in both Ẑ and Ŵ (Fig. 17), very similar to that observed in the synthetic model of an elevated 2-D conductor (cf. Fig. 5). This explanation seems to be particularly reasonable given the fact that in the original model the channeled current flows within the material with the conductivity σ ≈ 0.3 S/m and has the effective cross-sectional area of about ~ 10 4 m 2 (Fig. 4), which is equivalent to the current flowing in a long-length metallic ( σ ~ 10 6 -10 7 S/m) conductor with the cross-sectional area of just about several square centimeters.
On closer examination, it turned out that all 5 remaining non-causal impedances also reveal behavior very similar to that produced by the conceptual model in Fig. 4, which implies that the corresponding MT measurements be taken in a generally resistive environment with some elevated 2-D conductor enclosing a huge portion of local telluric currents. With that, 4 out of these five MT records were tagged by the field personnel as those acquired in a close proximity (< 300 m) to high-voltage overhead power lines, which in most cases are equipped with thick steel/aluminum ground-wire cables used for lightning protection. As for the last record from the "red" list of Table 2, we have no a priori information about the presence of power lines or other overhead metallic constructions in the survey area, but the latter is known to be relatively flat, so a reasonable geological explanation for an elevated conducting body could hardly exist anyway. The obtained results show that application of the above search algorithm to the whole onshore MT/ BBMT database of Nord-West Ltd., which includes thousands of records taken in complex geological conditions and extremely rugged topography (e.g., those acquired in Sub-Andean fold belt-see Palshin et al. 2020) has not yield a single impedance tensor, which could be confidently regarded as being non-causal due to geological reasons. This allows concluding that, though the potential existence of non-causal impedances (and hence magnetic tensors, tippers and any other transfer functions formally treating the horizontal magnetic field as an input-see Part I) on the land surface should not be excluded completely, the actual chance to encounter such data in practice is negligible.

Seafloor MT exploration
Offshore magnetotellurics, with hardly more than several hundreds of records annually acquired worldwide, is less popular than its onshore equivalent for the following reasons. The shallow-water (< 100 m) modification of the method, being highly vulnerable to the wave-motion noise and other sources of distortion (Constable et al. 1998), still remains rather undeveloped and only a few successful projects were reported recently (e.g., Geraskin et al. 2010;Ueda et al. 2014, Epishkin et al. 2018. The deep-sea MT exploration technique is equipped with state-of-the-art instrumentation, but the latter is very expensive and is deployed only from large sea vessels (Constable 2013). Moreover, though the receiving units are usually left underwater for many days of continuous operation, the effective bandwidth of the resulting data rarely exceeds 2-3 decades of period (starting from about 10 2 s) due to extreme attenuation of the high-frequency part of MT spectrum in saline water.
Among other things, the described specifics of marine magnetotellurics make the task of obtaining representative DR violation statistics in seafloor data rather challenging. However, hoping to reveal at least the main features of the offshore MT observations, we applied a similar search algorithm as used above to the integrated database of projects that Earthquake Research Institute (ERI) of the University of Tokyo and/or Japan Agency for Marine-Earth Science and Technology (JAMSTEC) participated in. This database includes: Normal Oceanic Mantle (NOMAN) project in northwestern Pacific (Baba et al. 2013(Baba et al. , 2017a, Stagnant Slab Project (SSP) in Philippine Sea (Seama et al. 2007;Baba et al. 2010), Mantle Electromagnetic and Tomography (MELT) experiment on East Pacific Rise (Evans et al. 1999;Baba et al.2006a, b), Tomographic Investigation by seafloor ARray Experiment for the Society hotspot (TIARES) project in French Polynesia (Nolasco et al. 1998;Tada et al. 2016), as well as marine data acquired in the vicinity of the Sado ridge (Toh et al. 2006), Mariana trough (Matsuno et al. 2010), and Tristan-da-Cunha hotspot (Baba et al. 2017b). The obtained results are given in Table 3. As to be anticipated, the seafloor MT measurements generally do not show any anomalous data far enough from coasts and other sharp bathymetric variations (e.g., NOMAN project), but may have plenty of such data when the observations Fig. 16 The DR application to the Q xy transfer function observed at site 28-002 in Yakutia, Russia are made in close proximity of continents, islands or seamounts (e.g., TIARES project). With that, vast majority of the revealed anomalous impedances were found to be non-causal. This fact suggests that the DR violations encountered in long-period offshore MT data are largely associated with irregular behavior of magnetic field due to intensive telluric currents flowing along underwater slopes.
By way of example we consider a seafloor MT record from the site Tris-11, acquired 3578 m underwater in the course of the cooperative German-Japanese academic project near the Tristan da Cunha archipelago (Baba et al. 2017b), just about ~ 25 km NE of the main island. The observed impedance components show prominently non-causal behavior for all true azimuth angles of the corresponding electric field except those about ~ 33° (and, consequently, ~ 213°), which could be easily recognized as the principal direction from the survey site to the nearest shore line (Fig. 18). This fact suggests that the same mechanism of "compensating zeros" as that securing the MP properties for the H-polarized (transverse) component of a non-causal tensor Ẑ in 2-D models (Fig. 8), applies to the H-polarized (radial) component of Ẑ in axially symmetric 3-D models as well.
The DR violation in off-diagonal components of Ẑ obtained at site Tris-11 is shown in Fig. 19. Though only Z yx reveals strongly pronounced non-causal behavior, both of these functions belong to the class 3 of the DR validity (Table 1). As a result, the dispersion relations should not be directly applied for quality assessment of Z xy as well, since it is difficult to tell whether the revealed minor DR violations in it are associated with actual data noise or merely with a slight mismatch of the corresponding lower half-plane poles and zeros (see Fig. 7). In accordance with Table 1, for this task one should find some other transfer function or a set of functions from the given response space, which would effectively carry the same information but be causal. It appears that for the offshore observations the most convenient choice of the kind is using the admittance tensor Ŷ = 1/Ẑ . Indeed, as follows from the above modeling and the present statistical review, the long-period MT measurements at the seafloor are mostly characterized by regular (MP) behavior of the horizontal electric field, which implies the causality of Ŷ . Figure 20 shows the off-diagonal components of the admittance tensor Ŷ representing the same MT data acquired at site Tris-11. Now, Ŷ is causal and Y yx is minimum phase (class 1), which allows to directly attribute any significant DR violation observed in Y yx to noise, source effect or other data quality issues. As for the Y xy component, with its phase rolling out of quadrant upwards (or, equivalently, rolling into the relevant quadrant from below) and DR-II violated by way more than ~ 180°, this function is recognized as one belonging to the class 2 of the DR validity with at least two non-MP zeros ( n ≥ 2 ). Consequently, correct application of the DR technique to such data can be realized either using the modified DR-II with several additional variables or, simply, the DR-I (Fig. 20, right).

Discussion and conclusions
Correct interpretation of the PROQ phenomenon, which used to be a daunting challenge at the times of simplistic 1-D and 2-D inversion codes (Jones et al. 1993), nowadays does not present a serious problem. Indeed, the electric field reversals responsible for non-MP behavior of the impedance components widely  (Table 2) are now treated as a useful constraint rather than a major impediment for MT interpretation (Ichihara et al. 2013), and can be effectively reconstructed by modern algorithms of 3-D inversion (e.g., Ichihara et al. 2014;Bologna et al. 2017;Pina-Varas and Dentith 2018). In the similar way, the magnetic field reversals responsible for the breakdown of the impedance and tipper causality in offshore MT data (Table 3) are readily accounted for by introducing detailed seafloor topography (Constable et al. 2009;Key and Constable 2011) and thus is no longer regarded as an obstacle-on the contrary, the distorted magnetic field is even reported to enhance the MT sensitivity to some geoelectric structures (Worzewski et al. 2012).

Fig. 19
The DR validity in a non-causal impedance tensor observed at the seafloor site Tris-11 (Tristan da Cunha) At the same time, the issues of appropriate processing and quality assessment of the MT data with violated dispersion relations still have some blind spots and open questions. In our research, we make significant progress on the subject by showing how to recognize the non-MP (class 2) and non-causal (class 3) behavior of a given MT response, and proposing three methods of the DR technique application to the revealed non-MP transfer functions: • employ the DR-I instead of DR-II, • employ the modified DR-II with n unknowns, or • try rotating the whole tensor unless all its significant components become MP. But which method should be chosen in practice? Since the DR-I is directly applicable to any component of any causal tensor, the first approach is clearly the most universal and is already recognized by some MT practitioners (e.g., Hacioglu et al. 2020). However, the DR-I application also implies linear vertical scaling of the functions under study. Though it might be advantageous for the MT transfer functions oscillating around zero value (such as the tipper or diagonal impedance components), linear scaling appears to be less favorable for the apparent resistivity curves since it effectively flattens the low-background anomalies while overemphasizing the high-background ones (see Fig. 16). In this context, we should expect that effective DR-I application to the MT data characterized by large range of values may require additional normalization procedures or other actions. The modified DR-II employs conventional log-scale representation of MT data, but its application also requires making use of n unknown variables (Eq. 7 in Part I). Though not a problem for n = 1 (Fig. 16), it might become an issue for n ≥ 2 (Fig. 20). The third method, based on the rotation procedure, is very straightforward but inapplicable to "strongly anomalous" MT tensors characterized by the presence of non-MP components at any rotation angle-such as the ones observed at the site B over the SL-shaped model (Additional file 1) or site 28-002 of the Chara-Tokko MT project (Additional file 2). Furthermore, the mere fact of using different coordinate axes for data processing, post-processing and inversion may lead to unnecessary quality losses. From this standpoint it seems advisable to keep the post-processing operations in accordance with the subsequent inversion routine, both in terms of the employed coordinate axes and data representation formats.
Another open issue of the DR application we see is the problem of quantifying errors. When employed at the post-processing step, the dispersion relations may effectively help to find less biased estimates (Fig. 5 in Part I), or even fill the gaps with no data at all (Fig. 8 in Part I). Apparently, this would improve the overall quality of the estimated transfer function, but for how much? Should the weighting matrix for the subsequent inversion routine be based upon the initial data error estimates at the given frequencies, or some special approach to the data weighting need to be introduced? Similar questions arise when the DR are applied for the quality assessment of the processed MT data. How to quantify the misfit between an actual MT curve and the corresponding DR value? Which norm should be chosen for the data represented as real and imaginary parts of a transfer function? And which one to use for a function given by phase and logarithmic amplitude, if the latter approaches infinite values (either negative or positive)? We suppose that finding correct answers to the raised questions would require further research and hands-on experience. However, some useful ideas on the subject could already be found in the papers discussing the general issue of choosing optimal misfit space for the inversion procedure (Miensopust et al. 2013;Wheelock et al. 2015;Miensopust 2017).
The main conclusion of our research is that it is legit and practically possible to apply the DR technique for assessing quality of virtually any acquired broadband MT data. If a function under study is causal, then any discrepancy between its imaginary part and the corresponding DR-I curve, as well as any discrepancy (other than related to the PROQ phenomena) between its phase and corresponding DR-II curve comes from poor data quality. If an MT tensor under study is non-causal, it can be equivalently represented via one or several causal tensors from the same response space.
The other principal conclusion refers to the old and heretofore unresolved dispute over the general existence of non-causal impedance tensors in real MT observations (see Berdichevsky and Dmitriev 2008). We show that these tensors do exist, but only at the sites with "abnormal" (non-MP) behavior of horizontal magnetic field characterized by one or more phase flips, such as those regularly observed off shore (Constable et al. 2009). As for the onshore measurements, though we succeeded to construct a couple of synthetic situations (one for a class of 2-D models with an elevated conductor, and another one-for a flat 3-D model) capable of producing a noncausal impedance function, none of them seems to be geologically reasonable. Furthermore, the subsequent examination of 60,000 + MT/BBMT onshore records yielded only a handful of non-causal tensors, with most, if not all, of them being non-causal due to the influence of nearby artificial conductors rather than some exotic geological conditions. This result is of high practical importance since it implies that the actual chance to encounter non-causal impedances, tippers and magnetic tensors in terrestrial MT exploration is vanishingly small. Given that the tipper function is especially vulnerable to noise and the source-effect distortions (Jones and Spratt 2002;Utada 2018), the relevant DR-I application for quality assessment of the tipper components proposed in (Marcuello et al. 2005;Shimizu et al. 2011) and justified by our research should be particularly useful.

Additional file 1. Synthetic data pack.
Additional file 2. Field data pack.