Magnetotelluric sampling theorem

We consider a general case of a magnetotelluric (MT) study to reveal three‑dimensional (3D) distribution of the elec‑ trical conductivity within the Earth based on measurements of electromagnetic (EM) fields by a two‑dimensional (2D) array. Such an MT array observation can be regarded as a spatially discrete sampling of the MT responses (imped‑ ances), and each observation site can be regarded as a sampling point. This means that MT array measurements must follow the Nyquist–Shannon sampling theorem. This paper discusses how the sampling theorem is applied to MT array studies and what kind of consideration is required in the application on the basis of synthetic model calculations, with special attention to spatial resolutions. With an aid of the EM scattering theory and the sampling theorem, we can show that an observation array resolves some features of the MT impedance but does not others. We call the resolvable and unresolvable features the MT signal and noise, respectively. This study introduces the spa‑ tial Fourier transform of array MT data (impedances) which helps us investigating sampling effects of lateral hetero‑ geneity from a different angle (in the wavenumber domain). Shallow heterogeneities cause a sharp spatial change of impedance elements near structural boundaries. High wavenumber Fourier components are required to describe such a feature, which means the site spacing must be sufficiently short to be able to resolve such features. Other‑ wise, a set of array MT data will suffer from aliasing, which is one of the typical causes of MT distortion (MT geologic noise). Conversely, a signal due to a deep‑seated conductivity anomaly will have more reduced amplitude at higher wavenumbers, which means focused imaging of such an anomaly is generally difficult. Finally, it is suggested to prop‑ erly consider the sampling theorem in an observation array design, so as to have best performance in resolving MT signals.


Introduction
Every geophysical exploration method usually conducts a sampling of physical quantities in the space domain.Sampled quantities are then translated into spatial distribution of a physical property in the Earth.For example, a geomagnetic survey samples magnetic field components in the space domain which are then inverted to a distribution of magnetization in the crust.Gravity measurements are inverted to density distribution in a similar manner.In these potential methods, the general Nyquist-Shannon sampling theorem is applied to argue the spatial resolution, the Nyquist wavenumber, etc.This paper concerns how the sampling theorem is applied to the practice of the magnetotelluric (MT) method (e.g., Cagniard 1953;Chave and Jones 2012).
The MT method is one of the geophysical methods to explore the Earth's interior based on measurements of time variations of electromagnetic (EM) fields on the Earth's surface or at seafloor.The MT response function (also called the MT impedance) is estimated at each observation site from time series of observed EM field variations as a function of frequency (e.g., Chave and Jones 2012).In the original derivation by Cagniard (1953), it was clearly demonstrated that the Earth's electrical structure can be estimated using the frequency dependence of the MT impedance obtained at an observation site, if the Earth is horizontally stratified, or in other words, one dimensional (1D).Because 1D Earth exists only in concept, an array of MT observations is necessary in general to reveal lateral variations in the MT impedance, so as to explore laterally heterogeneous Earth structures.
Time-changing electromagnetic fields follow Maxwell equations, which are written in the frequency domain as, where r =   x y z   , ω, µ 0 , and σ denote an arbitrary loca- tion below the surface ( z > 0 ), the angular frequency of time variations, the magnetic permeability of the vacuum, and the electrical conductivity, respectively.Note that the displace current is ignored in Eq. ( 2).
Due to the linearity of Maxwell equations and boundary conditions, the horizontal components of electric and magnetic fields measured at a certain location r h = x y at the surface of the Earth ( z = 0 ) are supposed to have a linear relationship with each other, where E h and H h are the horizontal electric and magnetic fields, respectively, both consisting of two horizontal components.Z is a complex-valued 2 × 2 tensor called the MT impedance.Now, a measurement at each observation site in the MT array is regarded as a sampling point of the MT (1) ∇ × E(r, ω) = −iωµ 0 H(r, ω), (2) ∇ × H(r, ω) = σ (r)E(r, ω), (3) E h (r h , ω) = Z(r h , ω)H h (r h , ω), impedance, and therefore, it is suggested that the observation array has to be properly designed on the basis of the general sampling theorem.Further, additional consideration is needed, because the spatial scale and time scale (frequency) of EM fields are related as will be discussed in later sections.Note that argument so far has been made by simply assuming all measured components completely satisfy Eq. ( 3), which means noise is absent.This does not happen in reality: any measurement includes noise at least to some extent.
One of the main concerns of this paper is how to deal with signals and noises in the MT method according to the sampling theorem.However, we need more clarification of the terminology of noise, because there are several different kinds of noise in MT practice including those in field observation, in data analysis and in inversion.One is the so-called measurement noise during observation, which is a contamination to original time-series data of electric and magnetic field variations.It could be caused by the measuring amplifier itself, by electrodes, by industries and railways, etc. (e.g., Ferguson 2012).In the inversion stage, inaccuracies due to physical approximations, truncations, or limitation of grid discretization could also be referred to as modeling noise.Among many different types of noise, this paper pays attention to the so-called geologic noise and corresponding (geologic) signal.
In this context, we first need a clear definition of 'signal' and 'noise' .Bahr (1991) tried to classify geologic noises in MT data into a number of types, but his definition of noise assumes the regional structure to be two dimensional (2D).In this paper, we first propose a definition of both geologic signal and noise on the basis of the sampling theorem, which is applicable to a general case of 3D heterogeneous Earth.

Signal and noise in MT and the sampling theorem
Here, we use the electromagnetic (EM) scattering theory (e.g., Hohmann, 1975;Chave & Smith 1994) to define signal and noise.The electric and magnetic fields, E and H, at an arbitrary position r in heterogeneous Earth can be formally expressed by integral equations as, (4) where E 0 and H 0 are the primary EM fields for a one- dimensional (1D) Earth σ 0 (z) which satisfy a condition of plane-wave approximation due to the incidence of an external EM field.G E and G H are the electric and mag- netic Green's tensor (e.g., Avdeev et al. 1997) and V k is the volume of each conductivity heterogeneity (anomaly) having conductivity contrast with respect to σ 0 (z) as For simplicity, σ is assumed to be constant in each vol- ume, V j .The second terms of the right-hand side of Eqs. ( 4) and ( 5) are called the secondary (scattered) fields that are laterally heterogeneous.
Although Eqs. ( 4) and ( 5) hold for arbitrary 'background' 1D profile σ 0 (z) , we take a regional mean 1D pro- file (Rung-Arunwan et al. 2017) that minimizes the lateral variance of δσ defined by Eq. ( 6).For a given MT obser- vation array, it is derived by an area average of the electrical conductivity at each depth slice either in a linear or in a logarithmic scaling, i.e., or where S denotes the area of the observation array.The averaged conductivity profile given either by Eq. ( 7) or by Eq. ( 8) is called the regional arithmetic or geometric mean 1D profile, respectively, which approximates σ 0 .
Let us consider a realistic observation where observed EM fields contain noise.In this case, horizontal components of the electric and magnetic fields measured at a location r h = x y at the Earth surface ( z = 0 ), denoted as E obs h and H obs h , can be formally written as, and ( 6) The first terms of right-hand side of Eqs. ( 9) and ( 10) are composed of horizontal components of EM fields given either by Eq. ( 4) or by Eq. ( 5).The second terms are those originated from a different (often very localized) source and, therefore, can be expressed neither by Eq. ( 4) nor by Eq. ( 5), which we call the measurement noise.
Next, we consider an MT array observation in which the impedance is estimated at each of discrete sampling points (observation sites).Here we assume a square array for simplicity with the length of one side L which we call the array size, and the number of sites on one side N (Fig. 1).Then, the site spacing is given by (10) which allows us to define the Nyquist wavenumber as It is well known that the resolution of wavenumber is given by 1/L .In this manner, we can directly apply the Nyquist-Shannon sampling theorem to the MT measurement.
The impedance, Z obs , at j th observation site r j = x j y j located at the Earth surface ( z j = 0 ) relates the measured EM fields in Eqs. ( 9) and (10) as, Here, δ represents contributions from measurement noises that do not follow the linear relation.Because measurement noises are involved in measured EM fields, the MT impedance has to be estimated statistically.
(12) We assume the effect of measurement noise in the impedance can be eliminated by a modern approach of statistical estimation (e.g., Gamble et al. 1979;Chave and Thomson 2004) and can be involved in the estimation error for impedance elements.Thus, Z obs can be regarded as an approximation of the impedance, Z , that relates horizontal EM fields without the effect of measurement noise as, In this paper, Z r j , ω is called the sampled impedance at r j .
On the basis of the Nyquist-Shannon sampling theorem, we can separate the scattered fields in Eqs. ( 4) and ( 5) into two parts: one whose spatial variations are resolvable by the given observation array and the other unresolvable.Here we call these resolvable EM field as MT signals.Significant effects from heterogeneities near the Earth surface with wavenumbers higher than the Nyquist and, therefore, not resolvable by the given observation array are defined as MT geologic noise, which are also called distortion.Ignoring measurement noise, the horizontal EM fields E h and H h sampled at a site r j = (x j , y j ) in the observation array can be formally written as, and where E 0 h are H 0 h are vectors consisting of two horizontal components of the primary fields in Eqs.(4) and (5) at the surface ( z = 0 ), e S and h S denote signals of the scattered horizontal electric and magnetic fields, and e N and h N denote their noises, respectively.
In the MT method, the EM fields without effects of distortion (geologic noise) are often called the regional (also undistorted) EM fields, which consists of the primary field and signal of the secondary fields, i.e., and where the superscript R in the left-hand side denote the regional fields.Then, we can define the regional (undistorted) impedance, Z R , as, ( 14) Revealing the electrical conductivity structure responsible for Z R is the main concern of the MT method, which obviously needs to remove the effects of geologic noises e N and h N contaminated in the observed EM fields (Eqs.15 and 16).

MT impedance elements in the spatial wavenumber domain
In this section, we examine features of EM fields and the MT impedance in the spatial wavenumber domain.The integral Eqs. ( 4) and ( 5) can be transformed to the spatial wavenumber domain at the Earth's surface (z = 0) as, where denotes the spatial Fourier transform and k x andk y are spatial wavenumbers in x-and y-directions, respectively.δ k E and δ k H denote the Fourier transform of the secondary electric and magnetic fields (the second terms of Eqs. 4 and 5) due to k-th conductivity heteroge- neity.For example, E x k x , k y , ω can be written as, By an analogy from a time-series analysis, we call the term with both wavenumbers equal to zero as the DC component and the rests as the AC components.These two equations indicate that the DC component provides the primary field and the AC components correspond to the secondary fields.Although it is possible to obtain the MT impedance in the spatial wavenumber domain from Eqs. ( 20) and (21) at the surface, it is not practical because all measurements have to be simultaneously carried out.
Therefore, we introduce the Fourier transform of the MT impedance to relate spatial distribution of the observed MT impedances to their spectra in the spatial wavenumber domain and to examine sampling effects (signal-noise relations) of the MT impedance.If all four (20 elements of the MT impedance at a certain frequency are known as a continuous function of position on the Earth's surface in a sufficiently large area, we can obtain the spatial Fourier transform of each element as, where Z ij is a (i, j)′ th element of the MT impedance.
If MT measurements are conducted at discrete sites distributed in a square array as shown in Fig. 1, we will have the discrete Fourier transform instead of Eq. ( 23) as, where l x and l y are integer between −N /2 and N /2 .The wavenumbers in x-and y-directions are given by l x /L and l y /L , respectively, both of which have discrete values between −k N and k N (Eq.12) with an increment of wave number resolution, 1/L .In the following sections, we call the dimensionless wave numbers l x and l y as the DLWNs.Z ij is the complex Fourier transform of Z ij , and we examine features of signals and noises based on wavenumber dependence of the amplitude spectrum,

Synthetic MT measurements in the spatial wavenumber domain assuming a checkerboard anomaly
In this section, we carry out numerical calculations for several synthetic models to examine features of the resulting MT impedances in the spatial wavenumber domain.The model space occupies 16 km × 16 km lat- erally and 300 km vertically, which was divided into 64×64×40 cubic cells.This model corresponds to the case of L = 16 km and N = 64 (see Fig. 1).EM induction equation was solved using a 3D forward code of WSIN-V3DMT (Siripunvaraporn, et al. 2005;Siripunvaraporn & Egbert 2009) at frequencies f = 100, 10, 1, 0.1, and 0.01 Hz, where f = ω/2π .Forward solutions provide us all components of electric and magnetic fields in response to two independent source fields and, therefore, 4 elements of the MT impedances at all grid points on the surface, as a function of frequency.Using these synthetic data at all grid points, one can calculate the spatial Fourier transform of the MT impedance from DC (23) ( l x = l y = 0 ) to the Nyquist DLWN ( l x = l y = 32 ), which we call the full spectrum hereafter.We examine features of synthetic MT data based on the spatial Fourier transform using checkerboard anomaly as a typical and simple example of 3D laterally heterogeneous Earth model which produces quasi-sinusoidal spatial variation of the impedance.We consider following 3 cases, in each of which the top of a checkboard with thickness of 0.1, 1.5, or 6.0 km is assumed at a depth of 0 (at the surface), 2 or 8 km, respectively (Fig. 2).Different checkerboard patterns are obtained by dividing each board into 2K × 2K (K = 1, 2, . . ., 16) square-shaped anomalies of alternating contrast of 1 order of magnitude in a uniform half-space of 0.01 S/m (Fig. 2).The symmetry of the model allows us to expect the impedance spectra with a fundamental mode at a spatial DLWN of l x = ±l y = K and its major higher harmonics also along l x = ±l y .Thus, we can examine major spectral structure only along an axis l x = l y for l x ≥ 0 .This simplicity of wavenumber domain characteristics is the main reason for using checkerboard models in the present synthetic study.

Features of full wavenumber spectra of MT impedances due to a checkerboard anomaly
The MT impedances were calculated from synthetic model with a checkerboard anomaly K = 2 .Figure 3 shows spatial distributions of the amplitude of impedance elements, Z xx and Z xy , at a frequency of 0.1 Hz for a checkerboard anomaly embedded at a depth of 2 km (those of Z yx and Z yy are not shown because of the sym- metry).Figure 4 shows spatial variation of the amplitude of Z xy along a line y = 2 (km) for a checkerboard anom- aly at each of three different depths (0, 2 or 8 km).We can see a sharp change of the impedance amplitude at the boundary of conductivity contrast when the checkerboard is exposed on the surface, but the spatial change becomes smaller and smoother with increasing depth.In the cases of the checkerboard depth of 2 km (Fig. 3), a clear checkerboard pattern can be seen in both elements, although the amplitude of Z xx is smaller than that of Z xy by nearly one order of magnitude or more.These impedances in the space domain are transformed to the spatial DLWN domain (Fig. 5).In the following, we mostly examine features of wavenumber spectra of the major element, which is denoted as Z xy (l x ) by making the fre- quency dependence implicit for simplicity.
Figure 6a shows the amplitude of full (wavenumber) spectra of Z xy (l x ) between zero ( l x = 0) and the Nyquist DLWN ( l x = N /2 ) for five different frequen- cies, and Fig. 6b shows its lower DLWN part ( l x ≤ 10 ) for more detailed inspection of features.We can clearly see spectral peaks at the fundamental DLWN ( l x = 2) and its higher harmonics up to l x = 6 , but the amplitude tends to decrease with increasing DLWN.This tendency reflects the gradual change in the impedance amplitude around the boundary between each of checker patterns (Figs. 3 and 4).
In Figs.7a and b, we compare the ratio of amplitudes at different DLWNs to DC amplitude ( Z xy (l x ) / Z xy (0) ).This normalization by DC amplitude allows us to see the wavenumber dependence more clearly, by canceling the f dependence of the impedance amplitude.For the fun- damental mode ( l x = 2 ), the value tends to be smaller at higher frequencies.This suggests that signals from conductivity anomalies attenuate more at higher frequencies through 2 km thick overburden with a cut-off frequency around 10 Hz. Figure 8a and b shows a similar plot for a case of the largest checker pattern ( K = 1 ) at the larg- est depth (8 km) examined in this study.We can see a rapid attenuation of signal ( l x = 1 ) amplitude at higher frequencies.In this case, the cut-off frequency is estimated to be around 0.1 − 1 Hz.Considering the depend- ence of the field penetration depth (induction scale length) given by Eq. (A4) in Appendix A on frequency, we can roughly estimate the ratio of the cut-off frequencies to be 1/16 , which is consistent with the estimations given above.
Further, we examine how full spatial wavenumber spectra of the MT impedances depend on the anomaly depth.Here, a checkerboard anomaly of K = 2 is consid- ered again, which is embedded at different depths, 0, 2, or 8 km.Resulting full spectra are shown in Figs.9a and  b.When the anomaly is exposed on the surface, spectral peaks can be clearly recognized at the fundamental DLWN ( l x = 2 ) and odd higher harmonic DLWNs up to the highest ( l x = 30 ) with gradually decreasing amplitude with increasing DLWN.All these higher harmonics can be ascribed to the sharp change in the MT impedance at the boundary of conductivity contrast in the space domain (Fig. 4).Peaks of even higher harmonics are also recognizable up to l x = 16 or so, which are supposed to reflect the asymmetry of impedance change due to a contrast in x-and y-directions (see Fig. 3).When the checkerboard is embedded at 2 km depth, the fundamental mode amplitude is similar to that of the previous case but higher harmonics can be traced only up to l x = 6 with much lower amplitude.When the checkerboard is buried as deep as 8 km, we can see only a peak at the fundamental DLWN with amplitude lower than the previous cases by more than 2 orders of magnitude.This indicates that such a checkerboard anomaly buried at 8 km depth is hardly resolvable unless the noise level is extremely low.
Figures 10a and b show full wavenumber spectra at a frequency of 0.1 Hz from 4 different checkerboard anomalies ( K = 1, 2, 4, and 8 ) embedded at a depth of 2 km.A peak at the fundamental DLWN can be recognized in each spectrum and the peak amplitude decreases with increasing K (smaller checker patterns).Peaks can be traced up to 2 or 3 higher harmonics for K = 1 and 2 , but not for other two cases.These features simply imply that anomalies with the larger size (the larger induction number) are better resolved.

Features of sampled wavenumber spectra of MT impedances due to a checkerboard anomaly
In the previous section, we examined major features of full wavenumber spectra, which is obtained from MT data sampled at every grid point.However, such a dense observation array in which each site locates exactly at each grid point is not available in actual practice of MT observations.Then, more realistic site distribution is considered for the synthetic sampled dataset, which consists of N × N sites distributed randomly (Rung-Arun- wan et al. 2017) as shown in Fig. 2. From 'sampled' MT measurements, we also calculate the discrete Fourier transform between DLWNs 0 ( l x = 0) and the Nyquist (l x = N /2 ) which we call the sampled spectrum hereafter.
First, we use a checkerboard K = 2 embedded at a depth of 2 km.Figures 11a and b show sampled spectra with different sampling intervals ( N = 32, 16, 8 and 4 ) together with full spectra (see Figs. 6a and b) at 0.1 Hz.The peak amplitudes at the fundamental DLWN are obtained rather accurately (close to that of the full spectrum), but it is biased upward a little for the smallest sampling number ( N = 4).At other DLWNs, spectral amplitudes are mostly estimated to be higher than the full spectrum amplitude.In other words, sampled spectra tend to more white for larger sampling intervals, which can be regarded as a typical feature of sampling effect.This feature explains the difficulty of imaging sharp structural boundaries of deep structures in MT inversion (Zhang et al. 2012).
Next, we sample the MT impedance for the same checkerboard ( K = 2) but exposed on the ground sur- face.As shown in Figs.7a and b., the full spectra of the Fig. 7 The same amplitude spectra as shown in Fig. 5 after normalized by the DC ( l x = 0 ) amplitude, showing larger attenuation of spectral amplitudes at higher frequencies calculated impedance with a sharp change at the exposed boundary consists of components from the lowest to highest DLWNs.Therefore, a strong aliasing effect from higher DLWN components appears in sampled spectra as shown in Figs.12a and b.Estimated amplitude values deviate further from that of the full spectrum for smaller sampling numbers ( N ).In particular, even the peak amplitude at the fundamental DLWN is estimated inaccurately.Inaccurate estimation of the peak amplitude eventually results in a false imaging of the underground structure.This example shows MT distortion can be caused by such a structure.Although this is not a typical example of MT distortion, we can learn from this result that careful treatment is required when there is an exposed structural boundary in the study area.
MT distortion is typically caused by anomalies which is smaller than the site spacing.Figures 13a and b show full and sampled ( N = 8 ) spectra of the MT impedances at 0.1 Hz from small checkerboard anomalies ( K = 8 and 16 ) exposed on the surface.In both cases of sampled spectra, the Nyquist DLWN ( l x = N /2 ) is smaller than the fundamental DLWNs ( l x = K ) indicated by arrows.As a result, the noise level of sampled spectra becomes higher than that of full spectra by 2 orders of magnitude.This is another simple example of aliasing effect (MT distortion).

Synthetic experiment of sampling effects of MT data with galvanic distortion
We further examine the sampling effects of distorted MT data.Here, we consider the case of purely galvanic distortion (see Appendix A).We use a simple parametric model of galvanic distortion (Rung-Arunwan et al. 2017), in which values of 4 independent parameters to describe the distortion operator C in Eq. ( 36) are given by random numbers with a normal distribution.Intensity of galvanic distortion is controlled by the standard deviation ( sd ) of the assumed normal distribution.We tried three different levels of sd = 0.02, 0.1, and 0.5 , which correspond to negligible, light and heavy distortions (Rung-Arunwan et al. 2017), respectively.The same sets of the sampled MT impedances as those in the previous sections are used as undistorted MT data Z R in Eq. (36) to generate a set of the distorted MT impedances CZ R .
Figures 14a and b compare features of full spectrum and those of sampled spectra obtained from the undistorted and distorted impedances at 0.1 Hz, using a checkerboard anomaly ( K = 2 ) embedded at a depth of 2 km.We can recognize spectral peaks at the fundamental DLWN in full and sampled spectra in cases without or with negligible (sd = 0.02) distortion.The spectral peak is masked by 'nearly white noise' in cases of heavier distortion, and the 'noise level' increases with increasing sd of galvanic distortion parameters.From this example, we recognize that 'MT geologic noise' is an appropriate terminology for galvanic distortion.
Next, we examine effects of sampling interval, using distorted MT data ( sd = 0.1 ) from the same checker- board anomaly ( K = 2, d = 2 km).We tested four cases of sampling number N = 4, 8, 16, and 32 , and result- ing wavenumber spectra are shown in Figs.15a and b together with the full spectrum for comparison.For the case of N = 32 , we can recognize a peak at the funda- mental DLWN with amplitude similar to the full spectrum.However, the peak is totally masked by noise for the cases of smaller values of N.
Here, we pay attention to the amplitude of DC components of the MT impedances at a frequency of 0.1 Hz obtained in previous sections.Interestingly, if impedances are not affected by galvanic distortion, the DC amplitudes of sampled spectra are nearly the same as that of the full spectra (see Fig. 6a), irrespective of the anomaly depth (Fig. 9a), the anomaly size (Fig. 10a), or the site spacing (Figs.11b and 12b).We can estimate similar amplitude for a case of weak or negligible distortion (Figs.14b and 15b) as well, but the estimated DC amplitude is significantly biased upward if the impedance is heavily distorted (Figs.14a and b).
Figures 16 a and b show the frequency dependence of sampled ( N = 8 ) amplitude spectra of undistorted (UD) impedance Z xy at a spatial DLWN l x = 0 and distorted ( sd = 0.02, 0.1, and 0.5 ) impedances Z xy at l x = 2 , for a checkerboard anomaly ( K = 2 ) embedded at a depth of 2 km.Obviously, the DC amplitudes are accurately estimated even for the case of severe distortion.The AC amplitudes, which mean the peak amplitudes at l x = 2 , are also accurate if distortion is not severe ( sd = 0.02 and 0.1 ).However, significant discrepancy appears at 100 Hz between the AC amplitudes of the undistorted and distorted impedances (undistorted amplitude is much smaller than distorted ones).This feature suggests that the MT signal from 2 km deep anomaly is too weak to be detected at this frequency (or higher) because of the strong attenuation already shown in Fig. 7. Overall, we have seen that estimates of the DC amplitude are more stable and accurate than AC amplitude estimates.When distortion is severe, DC amplitude estimates are slightly biased which may affect the reliability of the regional mean 1D profile.Rung-Arunwan et al., (2016;2017) pointed out that the rotational invariant MT impedance defined by ssq (Szarka and Menvielle 1997) of impedance elements, will give less biased estimate of the DC amplitude, instead of either of off-diagonal elements given by Eq. ( 24).A typical example of use of the regional mean 1D profile in 3D inversion can be found in Rung-Arunwan et al. ( 2022).( 27)

Discussion
It is generally difficult to handle aliased data.In case of a digital sampling in the time domain, this problem is solved by a use of an anti-aliasing filter.In the MT practice, there have been proposed several methods to correct the aliasing effects.For example, the electromagnetic array profiling (EMAP) is known as a field procedure that utilizes spatial filtering technique to avoid the spatial aliasing by a continuous sampling of the electric field (Torres-Verdin & Bostick 1992;Esparza & Gomez-Trevino 1996), which may be analogous to an anti-aliasing filter.However, its application is limited (Singer 1992) and not suitable for 3-D exploration.The Network-MT method by Uyeshima et al. (2001) reduces the influence of geologic noises by measuring the electric field with a large electrode separation, but its application is practically limited to areas where a network of metallic telephone lines is available (Hata et al. 2015).Bahr (1988) and Groom & Bailey (1989) proposed a scheme to correct the MT geologic noise assuming the distortion to be purely galvanic and the regional structure to be 2-D.Recently, full 3D inversion code with estimation of galvanic distortion tensor has been developed (Avdeeva et al. 2015;Usui 2015), but it is still a partial solution because site gains cannot be treated as a free parameter.The spatial spectral features shown in Figs. is when the MT impedances are heavily distorted.If distortion is inductive, the effect is expressed by a complex-valued tensor with frequency dependence (Utada & Munekane 2000), so that its removal becomes far more difficult than the case of purely galvanic distortion (Singer 1992).The presence of inductive distortion also violates the condition for the use of the MT phase tensor (Caldwell et al. 2004;Rung-Arunwan et al. 2022).In this context, a proper array design becomes an important issue as will be discussed below.
As shown in "Signal and noise in MT and the sampling theorem" section, configuration of an MT array (the array size L and site spacing L ) has to be designed by considering the necessary resolutions both in the space ( L ) and spatial wavenumber ( 1/L ) domains according to the sampling theorem.Moreover, additional consideration is needed in designing an MT array because the frequency and the spatial scale of EM induction are related.In the following example, we try to determine the frequency band of an MT array measurement so as to resolve all possible signals but not to be unnecessarily affected by geologic noises.
Let us denote the highest and lowest frequencies of the MT observation band as f max and f min , respectively.First, we consider the lowest frequency, f min , whose measure can be derived from a relation between the array size L and the maximum sounding depth required for resolving deep targets.The array size provides the maximum resolvable scale of the target structure to be L .Such a target structure embedded at a depth similar to its size is resolvable, if corresponding field penetration depth (induction scale length) at f min is sufficiently greater than the maximum resolvable scale.This condition can be written as, where 0 is the inductive scale length (skin depth) for σ 0 given by Eq. ( 35) in Appendix A and ε min is a number suf- ficiently smaller than unity.One can tune the value of f min by choosing a value of ε min .Equation ( 28) indicates that f min can be lower for a larger observation array and vice versa.It should also be noted that this equation casts an upper bound of f min of the observation band, instead of its lower bound.
In designing an MT array configuration, on the other hand, we must choose the observation frequency range so that any distorter ( d k < 2�L ) will not cause induc- tive distortion even at the highest frequency f max of the observation band (Singer 1992).This condition can be simply written as, (28) where M k and k are the induction number (Eq.34 in Appendix A) and the induction scale length (Eq.35 in Appendix A) of the corresponding distorter k , respec- tively, and ε max is another small number.Thus, we obtain a measure of f max with which one can expect a set of MT data free from being affected by inductive distortion as, Distortion is significant only when the conductivity contrast is large, and in such a case |δσ k | in Eq. ( 30) can be replaced by max[σ k , σ 0 ].
Figure 17 shows the relation between f max and L for typical values of max[σ k , σ 0 ] , assuming ε max 2 ∼ 0.1 .If we assume the same value for ε min 2 , the proper frequency band for the array study will be bounded by f max and f min as, Equation (31) suggests that the frequency range bounded by upper bounds of f max and f min may not be very wide, if considerations above were strictly applied.For example, the ratio f max /f min rarely exceeds 1000 even if the observation array consists of a large number of (29) (31) . Fig. 17 Relations between the site spacing, L , and the maximum frequency,f max , to avoid inclusion of inductive distortion for typical values of conductivity of 0.001, 0.01, and 0.1 S/m, estimated using Eq. ( 30) for a case of ε max 2 ∼ 0.1 sites, as the ratio σ 0 /max[σ k , σ 0 ] is at most 1.However, we should note that Eq. ( 28) gives the upper bound for f min rather than the lower bound, so that essentially no limitation is casted on f min .However, f min constrains the duration time t rec of data recording which is chosen to be sufficiently longer than the longest period 1/ f min ( t rec ≫ 1/ f min ).
In practice, MT observations are often made in a wide band (e.g., Ogawa et al. 1999), by letting the highest frequency be much higher than that given by Eq. ( 29), not by lowering the lowest frequency.This is a reasonable choice, because one can obtain more structural information on depths from wide band data without extending the recording time.However, one should note that such wide-band MT impedances may be affected by inductive distortion at higher frequencies.A careful examination about the type of distortion is indispensable.If the result obtained by 3D MT inversion suggests that the array design is not suitably optimized, one can improve it by performing additional measurements.Such a procedure is possible in MT study, because a set of impedances does not have to be obtained by simultaneous measurements.
Next, we show how the MT sampling theorem is applied to the synthetic modeling given in "Synthetic MT measurements in the spatial wavenumber domain assuming a checkerboard anomaly" section (L = 16 km,N = 64, σ 0 = 0.01S/m , |δσ k | ∼ 0.01 ).The array was designed to be applied to a local MT study.The vertical scale of the model space was chosen similar to the horizontal scale.The number of calculation grid is set to be sufficiently larger than practically possible number of sites, while it is set not too large to avoid too heavy computational burden for each model calculation.This model setup allows us to examine features of MT signals from a checkerboard anomaly of K = 1 ∼ 16 .Synthetic calcu- lation is made for discrete frequencies between f H and f L .The lowest frequency f L = 0.01 Hz was chosen so that the field penetration depths (induction scale length) to be comparable to the depth range ( 0 − 8 km) of the target structure (checkerboard anomaly).The site spacing ranges from 250 m ( N = 64 ) to 4 km ( N = 4 ), from which we can estimate f max (as a limit not to be affected by inductive distortion) roughly to be between 10 Hz and 4 × 10 −2 Hz, respectively, with ε max 2 ∼ 0.1 .Therefore, the highest frequency f H of the synthetic study is higher than this upper limit.Under this condition, we assumed in "Synthetic experiment of sampling effects of MT data with galvanic distortion" section that the synthetic MT data contain purely galvanic distortion.This means that we implicitly assumed the spatial scale of each distorter to be significantly shorter than its induction scale length (Eq.A3) even for the smallest value of f max so as all dis- tortions to be purely galvanic.In the present case, the typical scale of assumed distorters is supposed to be of order of 10 m or smaller.In this way, all distortions are set unresolvable with given values of the site spacing L .We finally optimized the model configuration so as to enable us to carry out all necessary inspections.It was possible in this synthetic study because we assumed distorters to be quite small, which we cannot control in an actual measurement.
In the case of seafloor MT array study to explore the oceanic mantle structure, a typical site spacing often exceeds 100 km, dealing with the highest frequency between 10 -2 and 10 -3 Hz.Considering the presence of significant bathymetric variations with scales smaller than the typical site spacing and the high conductivity of seawater, the condition of Eq. ( 30) is not satisfied in general and, therefore, observed impedances are supposed to be affected by inductive distortion.In fact, Baba et al. (2013) showed that the effect of bathymetric variations is expressed by a complex-valued distortion tensor especially at higher frequencies.Then, one might suspect that the use of the method proposed by Baba et al. (2010) is inappropriate, in which a model of the regional mean 1D profile is estimated by inverting the array-averaged (spatial DC) impedance.However in their method, impedances were averaged with iteratively updating correction of the bathymetric effects (distortion) at each site based on a forward calculation with known bathymetry.If this iteration converges, the impedance at each site is supposed to well approximate the regional (undistorted) impedance because distortion is accurately corrected.Thus, the estimation of 1D profile, as well as later estimation of 3D structures (Tada et al 2014), is justified in this case.

Conclusion
We considered how the Nyquist-Shannon sampling theorem is applied to MT studies, in particular, in a general case of the presence of distortion (geologic noise) over 3D heterogeneous Earth.Application of the sampling theorem clarified that signal and noise in MT data are defined as spatial features which are resolvable and unresolvable, respectively, by a given observation array.We applied the spatial Fourier transform of the MT impedance in order to examine behaviors of MT signals and noises in the spatial wavenumber domain, by using a set of synthetic MT impedances calculated from a checkerboard anomaly model.Major conclusions are as follows: (1) MT signals are features of impedances resolvable by the observation array, and those with smaller fundamental DLWNs (due to larger-scale anomalies) are better resolved.(2) MT geologic noise (distortion) is caused by a spatial aliasing of unresolvable (smaller than the site spacing of MT measurements) features of spatial variations of the MT impedance.(3) One can optimize the design of an MT array measurement in terms of array size, site spacing, and frequency range, by properly considering the sampling theorem and the basic physics of EM induction.

Appendix A. A brief review of galvanic distortion as a special case of MT geologic noise
An MT array measurement can be regarded as a sampling of impedance estimates at discrete points (observation sites).The given observation array can resolve spatial variations of the impedance due to laterally heterogeneous EM fields in a range of wavenumber between 1/L and the Nyquist (Eq.12).Because the geologic noise (distortion) is defined as the EM effects whose spatial scale is smaller than the Nyquist wavelength, it causes spatial aliasing (Torres-Verdin and Bostick 1992; Singer 1992) if noise amplitude is significant.Although the distortion is caused by some structures, it is neither possible nor interesting to estimate the responsible structures by inversion of the given dataset (a set of the observed impedances, Z ).Conversely, signifi- cant aliasing effects of distortion have to be properly eliminated like the case of Groom and Bailey (1989) to reveal the structures responsible for Z R by inversion within the allow- ance of estimation error.
In our definition, scattered EM fields whose lateral scale is smaller than the Nyquist wavelength ( 1/k N ) are regarded as MT geologic noise (distortion).Scattered electric fields may consist of both galvanic and inductive components, but scattered magnetic fields consist only of inductive (frequency dependent) component (Chave and Smith 1994).Therefore, the electric and magnetic noises in Eqs. ( 15) and ( 16) can be formally written as, where the superscripts G and I denote galvanic and inductive components, respectively.MT geologic noise due only to e G N is called galvanic distortion which is frequency independent.If both e I N and h I N persist, distortion is inductive and is frequency dependent.
We call near-surface small heterogeneities that cause distortion as distorters.It is not necessarily a lateral heterogeneity of conductivity, but topographic variations (Kaufl et al. 2018) can be a possible distorter.Utada and Munekane (2000) argued the significance of these noise (32) e N r j , ω = e G N r j + e I N r j , ω , (33) h N r j , ω = h I N r j , ω , components due to distorters of different spatial scales by simple order analyses as follows.Letting a typical scale of each distorter be denoted by d k , the induction number M k of each distorter can be defined as, where k is the induction scale length of k′ th distorter with conductivity contrast δσ k that can be written as, If M k ≪ 1 for a distorter, the contribution of inductive noise relative to the total electric field is roughly estimated as O(M k 2 ) , while it is O(M k ) relative to the total magnetic field.In other words, EM fields around such a distorter approach the static (galvanic) limit where time variation can be ignored.This means that the right-hand side of the Faraday's law diminishes with decreasing M k .In contrast, the galvanic term in Eq. (A1) little depends on M k but rather depends on the dis- tance from each distorter (more intense effect at shorter distance).If the induction number for all distorters contributing to the MT array measurement is significantly smaller than unity, inductive noise diminishes and distortion (noise) is purely galvanic.Thus, galvanic distortion can be regarded as a special case of MT geologic noise which is treated as purely galvanic at the frequency of interest.
Galvanic distortion can be formally expressed by a real-valued 2nd order tensor, denoted as C, operated to the regional impedance (Groom & Bailey 1989): i.e., the observed (distorted) impedance in Eq. ( 14) estimated at each site is related to the regional (undistorted) impedance as, D in the right-hand side of Eq. ( A5) is normalized so that where � • � F denotes the Frobenius norm (Bibby et al. 2005).The scalar normalization parameter g is called the site gain and the tensor D is called geometric distortion (Gomez-Trevino et al. 2013;Rung-arunwan et al. 2017), which is also called 'angular deviations of the induced telluric field' by Bahr (1991).In the MT method, g 2 is often called the static shift.Other parametrizations of galvanic (34) Z r j , ω = C r j Z R r j , ω = g r j D r j Z R r j , ω .
(38) g r j = �C r j � F , distortion have been proposed, but they are proved essentially equivalent (Smith, 1994).
Note that the distortion is galvanic and, therefore, C is real-valued and frequency independent, if and only if the induction number M k for every distorter is substan- tially smaller than unity for the entire frequency range of interest.Rung-Arunwan et al. (2016) proposed local or regional distortion indicator which will help judging whether such a condition is fulfilled for each site locally or in a given MT array dataset regionally.This condition is also necessary for the use of the MT phase tensor (Caldwell et al. 2004) as a distortion-free response function.

Fig. 1
Fig.1Configuration of a synthetic MT observation array, consisting of N × N square cells in each of which an MT observation site (black cross) is located.The cell size is L × L and the length of one side of the array is given by L(= N × �L) , which is called the array size

Fig. 2 A
Fig. 2 A plan view (right) of 3D checkerboard model in case of K = 2 used in the synthetic test.Checkerboard heterogeneity composed of 0.033 S/m (dark gray) and 0.0033 S/m (light gray) is embedded in the 0.01 S/m uniform half-space at a depth of 0, 2, or 8 km as shown in a side view (left) of the model.The right panel also shows an example of an observation array in case of N = 16 .Each sampled site location is denoted by a cross which is randomly located (Rung-arunwan et al. 2017) within each grid cell divided by dotted lines

Fig. 3 Fig. 4
Fig. 3 Synthetic MT datasets showing spatial distributions of Z xy (top) and Z xx (bottom) elements of the MT impedance at a frequency of 0.1 Hz due to a checkerboard anomaly of K = 2 buried at a depth of 2 km

Fig. 5 Fig. 6
Fig. 5 Amplitude of the spatial Fourier transform of Z xx (bottom) and Z xy (top) elements calculated from the synthetic dataset shown in Fig. 3

Fig. 8 Fig. 9 Fig. 10 Fig. 11 Fig. 12 Fig. 13 Fig. 14
Fig.8The amplitude of full spectra of Z xy normalized by the DC ( l x = 0 ) amplitude for the DLWNs l x (a) from 0 to the Nyquist and (b) from 0 to 10 at frequencies of 0.01, 0.1, 1, 10, and 100 Hz for a checkerboard anomaly K = 1 embedded at a depth of 8 km.The black arrow indicates the fundamental DLWN ( l x = 1) for the checkerboard anomaly

Fig. 15
Fig. 15 Differently sampled ( N = 4, 8, 16, and 32 ) amplitude spectra of distorted impedance element Z xy with sd = 0.1 at a frequency of 0.1 Hz for a checkerboard anomaly ( K = 2 ) embedded at a depth of 2 km, compared with the amplitude of a full spectrum for DLWNs l x (a) from 0 to the Nyquist and (b) from 0 to 10.The black arrow indicates the fundamental DLWN for the checkerboard anomaly

Fig. 16
Fig. 16Frequency dependence of sampled ( N = 8 ) amplitude spectra of (a) undistorted (UD) impedance element Z xy at a spatial DLWN l x = 0 and (b) distorted impedance element Z xy at l x = 2 , calculated for a checkerboard anomaly ( K = 2 ) embedded at a depth of 2 km.The distorted MT impedances were created using randomly generated distortion parameters with sd = 0.02, 0.1, and 0.5