Abstract
Electrical anisotropy has a significant impact on the observation data of the magnetotelluric (MT) method; therefore, it is necessary to develop forward and inverse methods in electrical anisotropic media. Based on the axis anisotropic electric field control equations, forming a large linear equation through staggered finite difference approximation, adding boundary conditions, and using the quasi-minimum residual method to solve the equation, this study obtained MT forward modeling results in axis anisotropic media. The correctness of the algorithm was verified by comparing it with the 2D quasi-analytic solution. By designing several sets of axis anisotropic 3D models, the characteristics of the apparent resistivity tensor and tipper were analyzed. The results indicated that the , and are sensitive to changes in resistivity in the X direction of the anomalous body, whereas the , and are sensitive to changes in resistivity in the Y direction. The apparent resistivity tensor and tipper are insensitive to changes in resistivity in the Z direction of the anomalous body. For exploration of anisotropic media, the apparent resistivity tensor and tipper of MT can identify the changes in resistivity in two horizontal axes directions and the boundaries of the anomalous body, which has the advantages for exploration.
1 Introduction
Anisotropy is commonly present in the crust and upper mantle. Fractures and rock bedding in certain specific directions, as well as the stacking combination of uniform thin layers with different properties, can cause electrical anisotropy in the lithospheric structure (; ). The practice and research of geophysical exploration have shown that electrical anisotropy has a significant impact on electromagnetic observation data and that directly using isotropic models to fit data containing electrical anisotropy can result in significant errors (; ). Therefore, to improve the accuracy of electromagnetic inversion results, including magnetotelluric methods, and the level of understanding of underground structures, it is necessary to develop electromagnetic data processing and inversion methods based on anisotropic models ().
Research on magnetotelluric anisotropy has been increasing (). Based on 2D numerical forward modeling techniques (; ), 3D forward modeling research results, especially the finite volume method (), the finite element method (; ; ; ; ; ) and the finite difference method (; ), continue to emerge. The finite element method has strong simulation ability for complex shapes and terrains, and it establishes variational equations through the Galerkin method. A weighted posterior error estimation method was constructed using the continuity condition of current density, which has been used to calculate the magnetotelluric (MT) response of arbitrary anisotropic media (). The finite difference method ensures that the distribution of the electromagnetic field satisfies the law of energy conservation and also simplifies the derivation of the equations.
Because of the complexity of numerical simulation for arbitrary anisotropy, this study considers the case of axis anisotropy. Axis anisotropy can be understood as the difference in conductivity in the three directions of the medium caused by factors such as mineral orientation. The conductivity tensor of arbitrary anisotropy can be obtained by three Euler rotations of the axis anisotropy. used direct discretization of Maxwell’s equations to achieve MT anisotropic forward modeling. In this study, the electric field control equation is solved to achieve MT axis anisotropic forward modeling. The algorithm is implemented by Fortran program. The response characteristics of the axis anisotropic target body are analyzed through three numerical examples, providing a basis for conducting MT forward and reverse modeling research with arbitrary anisotropy.
2 MT 3D forword modeling method
2.1 Finite difference method for calculating MT fields
For isotropic media, ignoring displacement current, the frequency domain control equation of the magnetotelluric method iswhere represents the Nabla operator, represents the electric field, represents an imaginary unit, represents the angular frequency, represents the vacuum magnetic permeability, represents the magnetic field, and represents the conductivity.
After organizing Equations 1, 2, the following electric field control equation is obtained:
For 3D axis anisotropic medial, the tensor conductivity is
By substituting Equation 4 into Equation 3 and organizing, the following equations can be obtained:
Using the staggered finite difference approximation Equations 5–7, the following equation is obtained after sorting ():where is a symmetric large sparse coefficient matrix, is a vector composed of three components of the electric field at the sampling point, and is a vector composed of boundary electric field components.
When a sufficiently thick air layer is added, the influence of anomalous bodies on the top boundary of the air layer can be ignored. The four lateral boundaries can be regarded as 2D geoelectric interfaces, solved using 2D MT anisotropic finite difference codes (). After introducing the boundary conditions, the three component values of the electric field in the partitioned space are obtained by solving Equation 8 using the quasi-minimum residual (QMR) method, with a preconditioner formed by an incomplete LU decomposition (). To accelerate the convergence of QMR iteration, divergence correction () is also applied to the solution of the electric field. The electric and magnetic field components of each measurement point on the surface are obtained through electric field interpolation.
2.2 Calculating tensor impedance and tipper
Calculate the magnetotelluric response using TE and TM polarization simulations, denoted as , , , and , , , , and . So the expressions for the magnetotelluric impedance component and tipper component are obtained:
The corresponding apparent resistivity tensor and apparent phase tensor are:where , . In the following text, the axial anisotropy response of MT is calculated by Equations 9–11.
2.3 Algorithm validation
To verify the correctness of the algorithm, we compared it with the quasi-analytic solution of , who established a 2D axis anisotropic upright fault model (Figure 1). The left side of the fault is an axis anisotropic block, with resistivity values of 10, 100, and 10 in the X, Y, and Z directions, respectively. The right side of the fault is an isotropic block, with resistivity values of 1,000 . A frequency of 0.1 Hz is used.
FIGURE 1
In TM mode, the comparison results of apparent resistivity and apparent phase and curve of relative residual error of the QMR iteration are shown in Figure 2. The 3D finite difference results and 2D quasi-analytic solution fit well, with only a slight error at the fault interface (Figures 2A–D), indicating that the calculation results of the forward program are correct. The QMR iteration converges stably to the given tolerance (Figure 2E).
FIGURE 2
3 3D axis anisotropic forward modeling case of MT
3.1 Response of 3D axis anisotropic low-resistance prism
The 3D prism model is shown in Figure 3. The top surface of the prism is buried at a depth of 100 m, and the prism has a size of 8 km × 8 km × 5 km. The background resistivity is 100 , and the mesh size is 28 × 28 × 28, which includes 7 air layers. The frequency of the MT used is 10 Hz.
FIGURE 3
To study the impact of the axis anisotropy of low-resistance bodies on MT forward modeling, we designed several sets of examples. We first fix the resistivity of the low-resistance body in the Y and Z directions ( and ), with resistivity in the X-direction of 10, 30, and 90 , respectively, and apply the algorithm proposed in this study for forward modeling. The resulting surface contour map is shown in Figure 4. Next, we fix the resistivity of the low-resistance body in the X and Z directions ( and ), with resistivity in the Y-direction of 10, 30, and 90 , respectively. These forward results are shown in Figure 5. Finally, we fix the resistivity of the low-resistance body in the X and Y directions ( and ), with the resistivity in the Z-direction of 10, 30, and 90 , respectively. The forward results for this example are shown in Figure 6.
FIGURE 4
FIGURE 5
FIGURE 6
The apparent resistivity and reflects the horizontal position of the low-resistance body, and the apparent resistivity , , and the tipper , clearly reflect the horizontal boundary of the low-resistance body (Figure 4). The amplitudes of and are much greater than those of and (Figures 4A–L). , , and are sensitive to changes in resistivity in the X direction () of the low-resistance body. When the resistivity of the low-resistance body is lower, the anomalies it produces are more obvious (Figures 4A–C, J–L, P–R). , , and are insensitive to changes in resistivity in the X direction of the low-resistance body (Figures 4D–I, M–O), and are less affected by the resistivity in the X-direction.
, , and are insensitive to changes in resistivity in the Y direction () of the low-resistance body, and are less affected by the resistivity in the Y direction (Figures 5A–C, J–L, P–R). , , and are sensitive to changes in resistivity in the Y-direction of the low-resistance body. As the resistivity of the low-resistance body becomes lower, the anomalies it produces become more obvious (Figures 5D–I, M–O).
The impacts of resistivity changes in the X and Y directions of the low-resistance body on response results are different.
The apparent resistivity tensor and tipper are not sensitive to changes in resistivity in the Z direction of the low resistivity body (Figures 6A–R), indicating that has a weak contribution to the surface response. For MT, in TE mode, the polarization direction of the electric field is mainly in the X direction, therefore it is sensitive to changes in resistivity in the X direction; in TM mode, the polarization direction of the electric field is mainly in the Y direction, therefore it is sensitive to changes in resistivity in the Y direction ().
A contour map of apparent resistivity based on the relationship between apparent resistivity and frequency is shown in Figure 7. The low-value anomaly areas of apparent resistivity and are not closed. As the detection depth of MT increases with decreasing frequency, the relatively low-value anomaly areas are equivalent to extending deeper underground (Figure 7).
FIGURE 7
3.2 Response of 3D axis anisotropic high-resistance prism
The 3D prism model is the same as 3.1, except that the resistivity is set to high-resistance. The frequency of the MT used is 10 Hz.
We first fix , , with of 200, 500, and 1,000 , respectively. The forward results are shown in Figure 8. Next, we fix , , with of 200, 500, and 1,000 , respectively. The forward results are shown in Figure 9.
FIGURE 8
FIGURE 9
The apparent resistivity and reflects the horizontal position of the high-resistance body, and the apparent resistivity , , and the tipper , reflect the horizontal boundary of the high-resistance body (Figure 8). , , and are sensitive to changes of of high-resistance body. When the resistivity of the high-resistance body is higher, the anomalies it produces are more obvious (Figures 8A–C, J–L, P–R). , , and are less affected by the (Figures 8D–I, M–O).
The apparent resistivity , , and the tipper are less affected by the of high-resistance body (Figures 9A–C, J–L, P–R). , , and are sensitive to changes of of the high-resistance body and are greatly affected by it (Figures 9D–I, M–O). The apparent resistivity tensor and tipper are not sensitive to changes of of the high-resistivity body.
3.3 Response of 3D complex axis anisotropic prism
The plan view of the 3D complex prism model is shown in Figure 10. The top surface of the model is buried at a depth of 100 m, the size of the prisms on both sides is 5 km × 22 km × 5 km, and the main axis resistivity in the X, Y, and Z direction is 10, 1,000, and 100 , respectively. The total length and width of the middle prism are both 12 km, and the height is 5 km. The main axis resistivity in the X, Y, and Z direction is 30, 500, and 100 , respectively. The background resistivity is 100 , and the frequency of MT is 10 Hz. The forward modeling results are shown in Figure 11.
FIGURE 10
FIGURE 11
Based on the previous analysis, according to Figure 11, the low-resistivity body in the X direction of the complex model caused the low value apparent resistivity anomaly zone of , and this anomaly zone reflects the boundary of the Y direction of the combination (Figure 11A). The high-resistivity anomaly zone of is caused by the high resistivity body in the Y direction of the composite, and this anomaly zone clearly reflects the boundary of the composite in the X direction (Figure 11B). The anomalous regions of and respectively reflect the boundaries of the combination in the X and Y directions, respectively (Figures 11C, D).
4 Conclusion
This study applies the finite difference method to achieve 3D MT forward modeling in axis anisotropic media and verifies the correctness of the algorithm by comparing it with 2D quasi-analytic solution. The examples show that , , and are sensitive to changes of of the anomalous body and are greatly affected by it, but insensitive to changes of . , , and are sensitive to changes of of the anomalous body and are greatly affected by it, but insensitive to changes of . The apparent resistivity tensor and tipper are not sensitive to changes of of the anomalous body. has a weak contribution to the surface response of MT. For the exploration of anisotropic media, the apparent resistivity tensor and tipper can identify the changes in resistivity in two horizontal axes directions and the boundaries of the anomalous body, which has the advantage for exploration. The study provides a basis for conducting the forward modeling and inversion research of MT with arbitrary anisotropy.
Statements
Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.
Author contributions
XL: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing–original draft, Writing–review and editing. Q-JS: Methodology, Validation, Writing–review and editing.
Funding
The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work is funded financially by Heilongjiang Province Basic Research Business Expenses for Universities Heilongjiang University Special Fund Project (Grant No. 2023-KYYWF-1494) and the Natural Science Foundation of Jiangxi Province (Grant No. 20212BAB213023).
Acknowledgments
Thank you to the reviewers for their valuable feedback on this article, Professor Siripunvaraporn, and others for their previous research work, and the editor for their enthusiastic assistance. Thank you to Wenxin Kong for his enthusiastic help. We would like to thank Editage (www.editage.cn) for English language editing.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1
CaoX. Y.YinC. C.ZhangB.HuangX.LiuY. H.CaiJ. (2018). A goal-oriented adaptive finite-element method for 3D MT anisotropic modeling with topography. Chin. J. Geophys.61 (6), 2618–2628. 10.6038/cjg2018L0068
2
GuoZ. Q.EgbertG.DongH.WeiW. (2020). Modular finite volume approach for 3D magnetotelluric modeling of the Earth medium with general anisotropy. Phys. Earth Planet. Interiors309, 106585. 10.1016/j.pepi.2020.106585
3
HanB.LiY. G.LiG. (2018). 3D forward modeling of magnetotelluric fields in general anisotropic media and its numerical implementation in Julia. Geophysics83 (4), F29–F40. 10.1190/GEO2017-0515.1
4
KongW. X. (2021). Research on the three-dimensional anisot ropic inversion of magnetotelluric data. Ph. D. thesis. Beijing: China University of Geosciences.
5
LiY. G. (2002). A finite-element algorithm for electromagnetic induction in two-dimensional anisotropic conductivity structures. Geophys. J. Int.148 (3), 389–401. 10.1046/j.1365-246x.2002.01570.x
6
LiuX.ZhengF. W. (2024). Axis anisotropic Occam’s 3D inversion of tensor CSAMT in data space. Appl. Geophys.10.1007/s11770-024-1076-9
7
LiuY.XuZ. H.LiY. G. (2018). Adaptive finite element modelling of three-dimensional magnetotelluric fields in general anisotropic media. J. Appl. Geophys.151, 113–124. 10.1016/j.jappgeo.2018.01.012
8
LiuY. H.YinC. C.CaiJ.HuangW.BenF.ZhangB.et al (2018). Review on research of electrical anisotropy in electromagnetic prospecting. Chin. J. Geophys.61 (8), 3468–3487. 10.6038/cjg2018L0004
9
PekJ.VernerT. (1997). Finite -difference modelling of magnetotelluric fields in two dimensional anisotropic media. Geophys. J. Int.128 (3), 505–521. 10.1111/j.1365-246X.1997.tb05314.x
10
PostmaG. W. (1955). Wave propagation in a stratified medium. Geophysics20 (4), 780–806. 10.1190/1.1438187
11
QinL. J.DingW. F.YangC. F. (2022). Magnetotelluric responses of an anisotropic 1-D earth with a layer of exponentially varying conductivity. Minerals12 (7), 915. 10.3390/min12070915
12
QinL. J.YangC. F.ChenK. (2013). Quasi analytic solution of 2-D magnetotelluric fields on an axially anisotropic infinite fault. Geophys. J. Int.192 (1), 67–74. 10.1093/gji/ggs018
13
SiripunvarapornaW.EgbertbG.LenburyY. (2002). Numerical accuracy of magnetotelluric modeling: a comparison of finite difference approximations. Earth Planets Space54, 721–725. 10.1186/BF03351724
14
SiripunvarapornaW.EgbertbG.LenburyY.UyeshimaM. (2005). Three-dimensional magnetotelluric inversion: data-space method. Phys. Earth Planet. Interiors150 (1–3), 3–14. 10.1016/j.pepi.2004.08.023
15
SmithJ. T. (1996). Conservative modeling of 3D electromagnetic fields. Part II. Biconjugate gradient solution and an accelerator. Geophysics61 (5), 1319–1324. 10.1190/1.1444055
16
WangT.WangK. P.TanH. D. (2017). Forward modeling and inversion of tensor CSAMT in 3D anisotropic media. Appl. Geophys.14 (04), 590–605. 10.1007/s11770-017-0644-7
17
WannamakerP. E. (2005). Anisotropy versus heterogeneity in continental solid earth electromagnetic studies: fundamental response characteristics and implications for physicochemical state. Surv. Geophys.26 (6), 733–765. 10.1007/s10712-005-1832-1
18
XiaoT. J.HuangX. Y.WangY. (2018). 3D MT modeling using the T-Ω method in general anisotropic media. J. Appl. Geophys.160, 171–182. 10.1016/j.jappgeo.2018.11.012
19
YeY. X.DuJ. M.LiuY.AiA. M.JiangF. Y. (2021). Three-dimensional magnetotelluric modeling in general anisotropic media using nodal-based unstructured finite element method. Comput and Geosciences148, 104686. 10.1016/j.cageo.2021.104686
20
YinC. C.WeideltP. (1999). Geoelectrical fields in a layered earth with arbitrary anisotropy. Geophysics64 (2), 426–434. 10.1190/1.1444547
21
YuG.XiaoQ. B.ZhaoG. Z.LiM. (2018). Three-dimensional magnetotelluric responses for arbitrary electrically anisotropic media and a practical application. Geophys. Prospect.66 (9), 1764–1783. 10.1111/1365-2478.12690
22
ZhouJ. J.HuX. Y.CaiH. Z. (2021). Three -dimensional finite-element analysis of magnetotelluric data using coulomb -gauged potentials in general anisotropic media. Pure Appl. Geophys.178 (11), 4561–4581. 10.1007/s00024-021-02882-0
Summary
Keywords
MT, axis anisotropy, response characteristic, finite difference method, forward-backward algorithm
Citation
Liu X and Sun Q-J (2024) A study of 3D axis anisotropic response of MT. Front. Earth Sci. 12:1454962. doi: 10.3389/feart.2024.1454962
Received
26 June 2024
Accepted
12 August 2024
Published
22 August 2024
Volume
12 - 2024
Edited by
Ying Liu, China University of Geosciences Wuhan, China
Reviewed by
Ningbo Bai, Henan Polytechnic University, China
Changhong Lin, China University of Geosciences, China
Updates
Copyright
© 2024 Liu and Sun.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Qi-Ji Sun, sunqiji422@126.com
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.