## 1 Introduction  
The magnetotelluric (MT) method is an important geophysical technique that utilizes natural alternating electromagnetic fields to investigate the electrical structure of subsurface media. In recent years, the effectiveness of the MT method has been well-demonstrated in resource exploration and geological disaster detection, especially in mineral imaging and geothermal exploration [1; 2; 3; 4; 5; 6]. With the growing complexity of study areas and the diversification of application environments, the MT method faces increasingly complex challenges, necessitating improved data acquisition and more appropriate inversion methods to achieve high-quality images [7; 8; 9; 10; 11]. For instance, in geothermal resource exploration, particularly in northern China [12; 13], the distribution of thermal storage is mainly controlled by regional faults. Fault zones are crucial for the transportation and distribution of hot water, as they fracture and saturate the deep Ordovician greywacke, creating channels for the upward conduction of Earth's deep thermal energy. Therefore, resolving fault zones is the most critical aspect of thermal exploration in this region. To identify faults beneath a cover, broadband MT surveys and long-period MT surveys can be conducted using high-resolution data to enhance the resolution of deep narrow anomalies [14; 15; 16]. However, MT inversion is inherently multi-resolution, and the inversion model often differs significantly from actual geological structures due to the limited resolution of the MT method and the influence of complex subsurface structures [17; 18]. Additionally, sparse low-frequency data points are insufficient to resolve deep fault structures [19; 20]. Moreover, different survey designs lead to varying interpretation results, affecting anomaly interpretation and geological assessment. To better understand the relationship between inverted anomalies and tectonics, and to improve the resolution and accuracy of target localization beneath cover, it is essential to design an optimal survey system and assess how line length, station spacing, observation bands, and model grid resolution influence the imaging of target tectonics. This target-oriented survey design can be achieved through forward modeling, allowing for the pre-determination of the subsurface structure and the evaluation of its resolution in MT inversion. Numerical simulation methods for MT forward modeling mainly include the Finite Element Method (FEM) [21; 22], the Spectral Method [23], and the Finite Difference Method (FDM) [20; 24]. The Finite Element Method is mostly used due to its high computational accuracy, adaptability to complex structures, and unconditional stability. Previous studies have shown that forward modeling plays a crucial role in optimizing survey design and improving the accuracy of data interpretation [25; 26; 27; 28; 29]. [PERSON] (2010) [30] considered the influence of the coastal effect on bottom magnetotelluric soundings in the shelf zone, using mathematical modeling of the magnetotelluric field into a two-dimensional medium to choose the most effective methods for calculating the fields near the coastal line. Based on magnetotelluric modeling, [PERSON] (2023) [31] applied the compressive sensing theory to MT acquisition and developed an efficient acquisition approach that uses an optimized station layout with sparse irregular patterns. This new acquisition approach can lead to significant savings in both the acquisition cost and operational time while maximizing the information obtained from a given number of stations. [PERSON] (2019) [32] applied genetic algorithms and a local methodology integrating the Gauss-Newton and Conjugate Gradient (GNCG) techniques to test one-dimensional inverse modeling of synthetic magnetotelluric data, leading to the development of a hybrid genetic algorithm for MT inversion of one- and two-dimensional synthetic data.