## 3 MT Method
### 3D MT Inversion  
The 3D inversion to transform the data into information for a preferred resistivity model was conducted using the WS3 DINVMT code [43]. The smoothest models were achieved by finding stationary points of an unconstrained functional \(U(\mathbf{m},\lambda)\):  
\[U(\mathbf{m},\lambda)=(\mathbf{m}-\mathbf{m}_{0})^{T}\mathbf{C}_{\mathbf{m}}^ {-1}(\mathbf{m}-\mathbf{m}_{0})+\lambda^{-1}\Big{\{}(\mathbf{d}-\mathbf{F}[ \mathbf{m}])^{T}\mathbf{C}_{\mathbf{d}}^{-1}(\mathbf{d}-\mathbf{F}[\mathbf{m }])-X^{*2}\Big{\}}. \tag{6}\]  
where \(\mathbf{m}\) is the resistivity model, \(\mathbf{m}_{0}\) is the prior model, \(\mathbf{C}_{\mathbf{m}}\) is the model covariance matrix, \(\mathbf{d}\) is the observed data, \(\mathbf{F}[\mathbf{m}]\) is the model response, \(\mathbf{C}_{\mathbf{d}}\) is the data covariance matrix, \(\lambda^{-1}\) is the Lagrange multiplier, and \(X^{*}\) is the misfit level.  
Figure 4: Examples of data processing for MT-01 and MT-03 data; (**a**) raw MT data and (**b**) data after cross-power selection and smoothing.  
The 3D MT inversion was performed for 26 MT datasets from 32 MT stations. We excluded six datasets because of their poor quality. In this 3D inversion, four response parameters were used: the real and imaginary parts of the off-diagonal impedances, \(Z_{xy}\), and \(Z_{yx}\). Several 3D MT inversion studies have used only the off-diagonal impedance tensor for the inversion process because the magnitudes of the diagonal components (\(Z_{xx}\) and \(Z_{yy}\)) were too small and the errors were too large [44; 45; 46; 47]. To reduce the computational loads, the data were resampled to include 10 period values out of 72-80 values from the 3.125 \(\times\) 10\({}^{-3}\) to 104 s period range.  
The impedance data were inverted using a 3D model comprised of meshes, where the smallest mesh dimension in the horizontal direction was 500 m \(\times\) 500 m. Outside the coverage area where the MT stations were distributed, the mesh dimensions were increased laterally by a factor of 1.2 to avoid the boundary effects of the modeling domain. The first vertical thickness of the ground was set to 100 m, followed by an increase in thickness in the  
Figure 5: Regional conductivity strike direction of the Limboto Lake—Pentadio area based on Groom–Bailey decomposition technique.  
vertical direction by a factor of 1.2. The total number of meshes in the \(x\), \(y\), and \(z\) directions was 48, 48, and 28, respectively. A homogeneous half-space initial model of 100 \(\Omega\)m was used in the inversion scheme.  
A final 3D resistivity model was obtained after 10 iterations and \(\sim\)120 h of computing time, with a root mean square (RMS) misfit of 2.685. Figure 6 depicts an example comparison between the observed and calculated data from MT-01 and MT-03 (for all stations provided in Supplementary Material 2. Figure 7 shows a comparison between the observed and calculated data in terms of the apparent resistivity and impedance phase maps, as well as their RMS values (for all periods provided in Supplementary Material 3). The highest and lowest periods involved in the inversion process were 3.125 \(\times\) 10\({}^{-3}\) and  
Figure 7: Comparison between the observed and calculated data for periods of 0.003 to 50 s. For the apparent resistivity maps, the lowest RMS value of 0.2 was yielded at a period of 0.0096 s, whereas the highest RMS value of 1.41 was yielded at 0.003 s. For the impedance phase maps, the lowest RMS value of 8.74 was yielded at 2.0408 s, whereas the highest RMS value of 14.88 was yielded at 50 s.