Necessity of Terrain Correction in Magnetotelluric Data Recorded from Garhwal Himalayan Region, India

The magnetotelluric (MT) method is one of the useful geophysical techniques to investigate deep crustal structures. However, in hilly terrains, e.g., the Garhwal Himalayan region, due to the highly undulating topography, MT responses are distorted. Such responses, if not corrected, may lead to the incorrect interpretation of geoelectric structures. In the present paper, we implemented terrain corrections in MT data recorded from the Garhwal Himalayan Corridor (GHC). We used AP3DMT, a 3D MT data modeling and inversion code written in the MATLAB environment. Terrain corrections in the MT impedance responses for 39 sites along the Roorkee–Gangotri profile in the period range of 0.01 s to 1000 s were first estimated using a synthetic model by recording the topography and locations of MT sites. Based on this study, we established the general character of the terrain and established where terrain corrections were necessary. The distortion introduced by topography was computed for each site using homogenous and heterogeneous models with actual topographic variations. Period-dependent, galvanic and inductive distortions were observed at different sites. We further applied terrain corrections to the real data recorded from the GHC. The corrected data were inverted, and the inverted model was compared with the corresponding inverted model obtained with uncorrected data. The modification in electrical resistivity features in the model obtained from the terrain-corrected response suggests the necessity of terrain correction in MT data recorded from the Himalayan region.


Introduction
The magnetotelluric (MT) method was first introduced by Tikhonov, Cagniard [1,2] to estimate the resistivity structure of the Earth's interior by analyzing the simultaneously measured components of natural time varying electric and magnetic fields on the surface of Earth. With the development of high-precision sensors, efficient data processing, modeling and inversion techniques, it is now effectively used for various applications related to subsurface resistivity structure delineation at shallow, intermediate and deeper depths [3][4][5]. It has successfully been used in the exploration of various Earth resources, e.g., geothermal, mineral, oil and gas resources, etc. [6][7][8]. In seismologically active zones, permanent MT monitoring stations are deployed to investigate the stability of the MT transfer function [9][10][11][12]. The changes in elastic properties due to high pressure and temperature are known to induce changes in the electric properties of rocks [13], which generates the precursor sign of earthquake preparation. Suppose the signal-to-noise ratio (SNR) is sufficiently good. In that case, the significant change in the conductivity of the subsurface will result in the shift of the MT transfer function [14]. By continuously monitoring MT responses, [15,16] predicted changes in the resistivity of subsurface structures (from 30 Ωm to 0.2 Ωm) after the Loma Prieta earthquake (M7. 1,1989) in California. The MT method is economic and efficient in comparison to the seismic method for deep crustal structure studies in difficult and highly undulating terrains such as in the Himalayan region [17,18]. However, the undulating topographic features modify the pattern of current flow, and therefore, affect both the electric and magnetic field components to different degrees. Thus, the MT response function becomes distorted when the observation site is located over or in the vicinity of undulating topography. Topographic distortion is estimated analytically [19] by assuming the Earth's surface is a single-valued, twice differentiable function f (x, y) with appropriate expansion for the surface magnetic field. The estimated corrections depend on the frequency and external radii of curvature of the Earth's surface at the observation point and, independent of the topographic profile, lead to accurate results. Several studies modeled the topographic effects on MT responses using 2-D and 3-D modeling [20][21][22][23]. The problem of distortion due to undulating topography and near-surface inhomogeneity is reviewed in [24,25]. A special-purpose 3-D FEM code was used for topographic correction and the interpretation of MT data recorded from a mountainous geothermal region [26]. The distortion tensor stripping-off technique was proposed by Larsen [27] and implemented in [28] for use in 2D topographic corrections of MT data. A comparison of computed topographic effects using two techniques: (i) the distortion tensor stripping-off technique and (ii) by correcting horizontal electric and magnetic field components was carried out in [22], and it was found that the two methods generally produce similar results. Over the last decade, several 3D FEM-based hybrid algorithms were developed that permit one to incorporate topography into the model [23,29,30] to avoid the distortion caused by the topography. The authors of [31,32] analyzed the influence of complex topography in MT-observed data from the northern Luzon Island, the Philippines and the Guangxi area, China. Using the 3D topography around the Hangai Mountains in Mongolia, topographic distortion was estimated [33] and elaborated in a general physical process within the model by analyzing current density variation and charge accumulation at the boundaries. The authors of [34] utilized the distortion tensor stripping-off approach to correct MT responses recorded from the Sikkim Himalayan region using 3-D modeling. They compared 1-D inverted models obtained using uncorrected and topographically corrected MT data recorded at selected sites in the Sikkim Himalayan region.
In the present paper, we applied the distortion tensor stripping-off technique to estimate topographic effects in the MT data recorded from the Garhwal Himalayan Corridor. For this purpose, we used MATLAB-based AP3DMT, a 3-D modeling and inversion code recently developed by our group [35]. The code was validated by estimating topographic effects over a standard trapezoidal model used in the literature [22]. We then used 3-D inversion on the topographic-effect-corrected MT data and compared the inverted model to that obtained from topographically uncorrected data. We noticed changes in the geometry of some of the geological features in the two models.

Methodology
We followed the distortion tensor stripping-off approach [27] for topographic correction, applied to the impedance tensor. In this approach, the distorted impedance tensor is linearly related to the undistorted impedance tensor through a distortion tensor, a twoby-two matrix (Equation (1)). If and are the topographically distorted and the flat surface undistorted impedance tensors, respectively, is the distortion tensor; the linear relationship is defined as: The distorted matrix elements of can be found by using two synthetic models: (i) with actual topography and (ii) the flat surface earth model. Both models must have the same subsurface resistivity structures. These models could be simple homogeneous half-space models with constant resistivity or may have a complex 3-D resistivity structure. After obtaining the distortion tensor matrix elements of through synthetic modeling, using Equation (1), the undistorted impedance tensor can be calculated from the recorded impedance response.

Validation of Code
To validate our AP3DMT code for 3-D response computations for an undulating topographic model, we used the trapezoidal hill model of [22] and replicated the results published in [22]. The trapezoidal hill model ( Figure 1) consists of 0.45 km-high hill that is 0.45 × 0.45 km wide at the top and 2 × 2 km wide at the bottom. An anomalous body of 0.388 × 0.388 km dimension is inserted into the hill from surface relief to a 2 km depth. The dimension of the model domain is 9.8 × 9.9 × 22.2 km, and we used 12 MT profiles consisting of 12 sites in each profile. Thus, a total of 144 sites data were used to study the response characteristics. For comparison, the grid, used in x-and y-directions, was the same as in [22], and in the z-direction, the grid spacing started from 50 m and it increased by a factor of 1.2. The elevation grid spacing used was, from surface to 450 m height, constantly 25 m. In addition, eight air gridlines were also included. The forward response computations were performed with topography as well as flat earth models, and the topographic corrections were computed using the method discussed in the previous section. Figure 2 shows the comparison of computations to that of published results [22] for the profile lines 6, 8 and 10 ( Figure 1 We further investigated the topographic effect over the entire model domain at 22 Hz. This was performed by interpolating the responses between two consecutive MT sites ( Figure  3). In this model, we used the topographic dyke (TD) of 10 Ω m resistivity placed in a 100 Ω m half-space. The low apparent resistivity values (less than 100 Ω m) are observed over the hill, and the region where conductive body is present and has high apparent resistivity values (more than 100 Ω-m) can be seen in the foothill regions and at the boundaries of the low resistivity structure. The TD responses show the variation in apparent resistivity between 3.6 to 240 Ω m and in the phase between 40° to 58°. This is in agreement with the known fact that electric fields are reduced on topographic hills and increased in valleys due to the focusing of current in opposite the sense in the two cases, and this is responsible for generating the galvanic effect [24]. The apparent resistivity values are high in valley troughs and low at the peaks. Therefore, the changes in apparent resistivity seen over the hill and valley regions are purely galvanic in nature. The distortion in TE-and TM-mode responses are consistent with the topographic features [33]. The topographic dyke (TD), corrected (TC) together with the flat-earth's dyke (FD) responses, is shown in Figure 3. Our next step was to investigate the effect of ramp, slope and distance from the slope. For this purpose, a topographic ramp model with different slope angles was used. Computed topographic effects at different points of the slope, hilltop and distant points from the slope are shown in Figure 4. MT responses show opposite effects in amplitude at sites located at the upside and downside of the slope, respectively. The topographic effects of the slope's angle at three frequencies (100 Hz, 22 Hz and 1 Hz) are shown in Figure 4b. The response shows that the effect is proportional to the slope's angle and that the topographic effect is negligible at the sites located more than 1 km from the slope. In the next section, we show how we applied topographic corrections to the data recorded from the GHC region.

Terrain Corrections in the Field MT Data from the Himalayan Region
Our group recorded broadband MT data along the Roorkee-Gangotri (RG) profile in the Garhwal Himalayan Corridor (GHC), Uttarakhand, India. The study area is located between 29.0° N to 31.5° N and 77° E to 80° E, and the RG profile passes across the major Himalayan thrusts in the GHC. The details of equipment and processing techniques used are given in earlier publications [36,37]. Impedance tensor responses estimated for 39 MT sites in the period range of 0.01 to 1000 s are used in the present work. Locations of the MT sites are shown in Figure 5. These sites are located at different elevations, starting from the Indian-Gangetic Plain (IGP) in the southernmost end to the Higher Himalayan (HH) crystalline region in the northern end of the profile. This profile passes through the highly undulating terrain and crosses four major concentric litho-tectonic domains, separated by three south verging major Himalayan thrusts. From south to north, these thrusts are the Main Frontal Thrust (MFT), Main Boundary Thrust (MBT) and Main Central Thrust (MCT). In the GHC, the MCT is defined as a zone bounded by the Munsiari Thrust in the south and the Vaikrita Thrust in the north, also referred as MCT-1 and MCT-2, respectively [38,39]. Here, we used the latter convention. Four litho-tectonic domains from south to north are the Indian-Gangetic Plain (IGP), the Sub-Himalaya (SH), Lesser Himalaya (LH) and Higher Himalaya (HH) domains. The IGP domain has an average elevation of 218 m in the study area. The IGP is separated from the Sub-Himalayan (SH) domain by the MFT in the north and the peninsular shield towards the south. The SH region begins to the north of the MFT with an average elevation of 600 m, and it rises abruptly above the IGP along the MFT. The Lesser Himalayan (LH) domain, which is located to the north of the MBT, has an average elevation of 2500 m, with an increasingly elevated profile towards the north. The LH domain is further divided into two zones: the southern zone, lying between the MBT and Sri Nagar Thrust (SNT), is the Outer Lesser Himalaya, and the northern zone, lying between the SNT and MCT-1, is the Inner Lesser Himalaya.  [40,41] The northernmost region of the profile is in the Higher Himalayan (HH) crystalline zone, containing most of the famous peaks of the mountain range, and it has an average elevation of 4500 m. In the main Himalayan arc, the dominant structural trend of the MFT, MBT and MCT are generally aligned in the northwest direction. The MT site's elevation along the RG profile varies between 218 m in the IGP to 3300 m in the MCT zone. In addition, this profile passes through the highly undulating Himalayan terrain and, therefore, the terrain corrections estimation is a logical intermediate step for this dataset before inversion. However, terrain correction has not been implemented in previous 3-D inversion studies [37]. In the following, we rexamined the entire data set and applied terrain corrections wherever needed and subsequently used 3-D inversion to the terrain corrected data set.
Terrain effects were estimated for all 39 MT sites of the RG profile using the impedance tensor stripping-off technique discussed in Section 2. We used two models: model-1, a homogeneous half-space model with resistivity of 100 Ω m and model-2, a heterogeneous, checkerboard model [35,42]. The heterogeneous model-2 consists of alternative blocks with resistivities of 1 Ω m and 1000 Ω m, respectively. The top of each block is located at a depth of 2 km from the flat earth surface (Figure 6). The dimensions of each block are 40 × 40 × 2.5 km in the x-, y-and z-directions, respectively, and the blocks are placed symmetrically in x-and y-coordinates at −80 to −40 km, −20 to 20 km and 40 to 80 km. Actual topographic elevations corresponding to sites in the RG profile were used. Data regarding elevations above mean sea level (MSL) and topographic variations were extracted from the Indian Geo-Platform of ISRO [43]. The model domain is 200 × 200 sq. km with 100 Ω m half-space resistivity used for modeling. Fine grid discretization was used near the MT site (within 2 km) for better resolution. Therefore, the x-and y-grid spacing was chosen in an increasing order from 20 m to 440 m in 4 × 4 sq. km area for each site zone. Within this region, a constant grid spacing of 50 m was used in the z-direction above the flat earth surface. The total height of the model depended on the site's elevation. The z-grid spacing started with 50 m and increased by a factor of 1.2 in depth below the flat earth's surface. The forward MT response was computed at 39 sites of the RG profile. The apparent resistivity and phase deviations from its values for the flat earth model were the topographic effects. For the homogeneous model, the changes in apparent resistivity beyond 100 Ω m and the changes in phase beyond 45° were purely due to the terrain effects. On the other hand, in the case of the heterogeneous model, the changes with respect to the flat-surface model response were due to the topographic effect. The estimated topographic corrections for both the models are presented in Figure 7 with regard to TE-and TMmode responses, respectively. It is observed that the nature of topographic effects is galvanic and both galvanic and inductive at different sites located in different litho-tectonic domains. The nature of the topographic effect is summarized for all sites in Figure 5a. The southern zone of the profile is located in the IGP, with an average topographic elevation of 200 m and no severe topographic undulation toward the southern end of the profile. However, the topography in the northern zone of the IGP is influenced by the MFT, the southernmost Himalayan thrust. Thus, the topographic effects observed in sites 06 are the results of undulation present in and around this site. To demonstrate the nature of topographic effects in the GHC region, we selected representative sites: 06, 14, 19, 31 and 40 located in the IGP, SH, Outer-LH, Inner-LH and HH regions, respectively, for detailed discussion in the following sections.  Site 06 is located at an altitude of 519 m, and the maximum variation in the elevation is along the N 43° W direction; the elevation varies from 507 m to 698 m, with the slope ranging between −35° and +13° within 1 km from this site. Topographic elevations along the NS and EW directions around this site are shown in Figure 8b. The site is located on a valley slope. Therefore, the topographic effect of both slope and valley contributed to the results. However, the dominant topographic effects in off diagonal components are mainly due to the valley and are shown in Figure 8d,e for the homogeneous and heterogeneous models, respectively. Topographic effects show different magnitudes of positive shifts in amplitudes of both (TE and TM) modes of polarization. In the literature, this is also referred to as static shift effect, which is mainly due to the accumulation of electrical charges at the interfaces of different electrical resistivity zones.    Figure 11b,c, respectively. The site is located near a river channel in the valley. Both galvanic and inductive topographic effects are observed at this site. However, the galvanic effect is dominant in the MT responses, as shown in Figure 11d. The amplitude of the two modes shows an opposite shift: upward in the TE-and downward in the TM-mode, respectively. The last selected site, site 40 ,is located in the Higher Himalayan region at an altitude of 2759 m. The elevation and slope vary from 2493 m to 3425 m and −57° to +35°, respectively, within 1 km distance from this site. The site is located on a hillslope. Elevation variations along the NS, EW and the maximum variation direction are shown in Figure  12b,c. Both galvanic and inductive effects can be seen in Figure 12d at this site. Thus, out of the total of 39 sites, we observed terrain effects in 25 sites. The nature of terrain effects is purely galvanic in 16 sites (blue circle, Figure 5), whereas terrain effects are both galvanic and inductive in 9 sites (red circle, Figure 5). Significant (>50%) terrain effects in amplitudes of and components are observed in 10 sites (6,14,16,24,34,35,37,40,44,52) and in 8 sites (15,19,20,32,34,36,37,49), respectively. Terrain effects in phase ( ) responses are observed in the lower period range (<10 s) and are relatively smaller in magnitude; they vary from 2 to 6 degrees at 11 sites (14,16,24,31,32,34,35,37,40,44

3D Inversion
We performed 3-D inversion of terrain effects for corrected and uncorrected full impedance tensors (Z) using our code, AP3DMT [35], and compared the two inverted models. The code was run on a HP Z620 workstation (with two E2643 3.30 GHz Intel Xeon processors and 32 GB RAM with eight cores). AP3DMT, a MATLAB-based code for the 3-D modeling and inversion of MT data, uses numerical library functions and subroutines efficiently [35]. The inversion algorithm AP3DMT is based on finding a meaningful Mdimensional model parameter vector, by fitting the -dimensional data vector, to an acceptable level in a stable manner. The program uses the minimization of a penalty functional, defined as , where is the forward mapping operator, Cd is the data covariance matrix, is the apriori model, Cm is the model covariance matrix or smoothing matrix and λ is the tradeoff parameter. The convergence of the iterative process is defined by the normalized Root Mean Square (nRMS) value, given as where N is the total number of data points; is the observed data; is the synthetic data and is the standard error in observed data value. The AP3DMT code allows the user to define (i) the model covariance matrix (Cm), (ii) the initial values of trade-off parameter (λ) and (iii) the decreasing factor for λ needed while updating the model parameters. Parameters defining smoothing can be chosen between 0 and 1. Inversion was performed by various test runs on the uncorrected dataset [37] and, therefore, we did not repeat the same exercise here; instead, we chose to set the inversion control parameters to the same as in [37] for both the uncorrected and corrected full impedance tensor data. The initial guess model is a 100 Ω m resistivity of a homogeneous half-space central domain of 200 × 200 sq. km. The model is discretized in x-and y-directions with a constant grid spacing of 4 km, and in the z-direction, the grid spacing starts with 50 m and is increased with a 1.2 factor in depth. Thus, the total number of grids in x-, y-and zdirections are 66 × 63 × 44, respectively. The error floor is set as 10% of | * | / in all the components. The isotropic smoothness is 0.3, and the starting regularization parameter (λ) is set to 10. In the case of the inversion of uncorrected data, the starting value of nRMS is 27.9, which converges to 2.006 in 81 inversion iterations. The inversion process was terminated at the regularization threshold criteria. For topographic corrected data, the initial nRMS is 30.97, which converges to the final nRMS value of 2.58 in 92 iterations.
In an attempt to investigate the reason for the increase in final nRMS values, we found some spikes and outliers in the data which were poorly fitted; this might have resulted in an increase in nRMS. In the next inversion run, we removed these spiky and outlier data points and ran 3-D inversion with an initial model as the final inverted model of previous run. This inversion converged to nRMS 2.25 in 32 iterations. We concluded that the inversion of uncorrected data partly fitted the error introduced in the response due to the terrain effect. Thus, the difference in nRMS in the two inversion runs was due to the fact that the inversion was over-fitting the response in the case of uncorrected responses. However, the model was better resolved in the case of the 3-D inversion of terrain-corrected data in spite of slightly higher nRMS. This point is further elaborated in the following detailed discussion on the two models. The two inverted models are plotted in the form of columns representing the depth slices and the ratio of cell resistivities in Figure 13a-c. This helps in understanding the regions of the same (ratio = 0) or different (ratio ≠ 0) resistivities in the two models. It can be seen that the cell resistivities ratio is zero in most of the region, except for the region where terrain effect is severe and capable of modifying the model resistivity, e.g., in the northern region of the profile. The high ratio value is observed at the NE side of the profile, where the topography is highly undulating. We further compared the profile sections of the two models, obtained from the terrainuncorrected and corrected responses, shown in Figure 14. A significant change in the two models is observed in the NE region of the profile. The resolution and geometry of the two geological features vary significantly in the two models. A conductive body ( 5 Ω ) beneath the Higher Himalayan zone is observed at a depth of about 3-6 km; it is denoted as the mid-crustal conductor, inferred as metamorphic fluids released due to continuous underthrusting of the Indian plate. This feature is common in the Himalayan region and was also observed by other researchers in other sections of the HH regions [44,45]. Another significant change is shown by the red dashed line (Figure 14), termed the detachment surface that has the flat-ramp-flat geometry of the Main Himalayan Thrust (MHT) in the GHC region. In the terrain-corrected model, this feature is clearer and more consistent. The flatramp-flat geometry of the MHT is consistent with other studies [46,47] carried out in the region.

Conclusions
We implemented terrain corrections in MT data recorded from the GHC region. We used the distortion tensor stripping-off technique implemented in AP3DMT, a MATLAB-based 3-D modeling and inversion code previously developed by our group. Our methodology was validated by replicating the results over standard models available in the literature. Subsequently, terrain corrections were applied in the MT data acquired along the RG profile in the period range of 0.01 to 1000 s. The distortion introduced by topography was computed for each site using homogenous and heterogeneous models with actual topography at the recording sites. Significant terrain effects were observed in amplitude as well as in phase responses, especially in the sites located in the LH and HH region. The models obtained from 3-D inversion of the terrain-corrected and uncorrected data were compared. This study suggests that terrain correction is necessary in MT data recorded from the GHC region, as terrain-corrected responses modify inverted model features. Code availability (software application or custom code): Research material and code can be shared with approval for research purposes.