- Email: [email protected]
2-D magnetotelluric modeling using finite element method incorporating unstructured quadrilateral elements
Journal of Applied Geophysics 139 (2017) 16–24
Contents lists available at ScienceDirect
Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo
2-D magnetotelluric modeling using finite element method incorporating unstructured quadrilateral elements Weerachai Sarakorn Department of Mathematics, Faculty of Science, Khon Kaen University, Khon Kaen, Thailand
A R T I C L E
I N F O
Article history: Received 25 May 2016 Received in revised form 10 February 2017 Accepted 13 February 2017 Available online 16 February 2017 JEL classification: 02.70.Dh 91.25. − r 2000 MSC: 65N30 65N50 74S05
A B S T R A C T In this research, the finite element (FE) method incorporating quadrilateral elements for solving 2-D MT modeling was presented. The finite element software was developed, employing a paving algorithm to generate the unstructured quadrilateral mesh. The accuracy, efficiency, reliability, and flexibility of our FE forward modeling are presented, compared and discussed. The numerical results indicate that our FE codes using an unstructured quadrilateral mesh provide good accuracy when the local mesh refinement is applied around sites and in the area of interest, with superior results when compared to other FE methods. The reliability of the developed codes was also confirmed when comparing both analytical solutions and COMMEMI2D model. Furthermore, our developed FE codes incorporating an unstructured quadrilateral mesh showed useful and powerful features such as handling irregular and complex subregions and providing local refinement of the mesh for a 2-D domain as closely as unstructured triangular mesh but it requires less number of elements in a mesh. © 2017 Elsevier B.V. All rights reserved.
Keywords: Magnetotelluric modeling Finite element method Unstructured quadrilateral mesh Paving algorithm
1. Introduction From the past to the present, two-dimensional magnetotelluric modeling was successfully solved by various numerical techniques such as the finite difference (FD) method (Rao and Babu, 2006) and finite element (FE) method (Reddy and Rankin, 1975; Wannamaker et al., 1986; Key and Weiss, 2006; Franke et al., 2007; Lee et al., 2009). The efficiency and accuracy of the FD method were proven for a simple 2-D model (Rao and Babu, 2006). With the limitation of rectangular grids, many complex subregions in the 2-D domain such as modeling topography, irregular anomalous bodies, and seafloor cannot be handled by its discretization. Therefore, the solution obtained by the FD method may be inaccurate and unrelialable. The finite element (FE) method is another numerical technique that is often used for 2-D MT forward modeling. Its efficiency and accuracy have been proven (Reddy and Rankin, 1975; Wannamaker et al., 1986; Key and Weiss, 2006; Franke et al., 2007; Lee et al., 2009). Unlike the FD method, the FE method is more flexible on discretization. Various mesh schemes such as structured triangular meshes (Wannamaker et al., 1986), unstructured triangular meshes with
http://dx.doi.org/10.1016/j.jappgeo.2017.02.005 0009-2541/© 2017 Elsevier B.V. All rights reserved.
adaptive refinement (Franke et al., 2007; Key and Weiss, 2006) and structured quadrilateral meshes (Reddy and Rankin, 1975; Lee et al., 2009) can be used in the FE method. The structured triangular and quadrilateral mesh algorithms are simple and can handle domains with complex subregions, but the local refinement mesh feature cannot be implemented efficiently. In contrast, the unstructured triangular mesh algorithm can overcome this problem. However, it needs to incorporate an adaptive scheme to improve its accuracy and it requires a large number of triangular elements in complex subregions (Key and Weiss, 2006; Franke et al., 2007). Presently, implementations of unstructured quadrilateral meshes already in the FE method are popular as discussed by Blacker and Stephenson (1991); Bern (1997); Owen (1998); Sarrate and Huerta (2000); Remacle et al. (2012). The advantages of unstructured quadrilateral mesh algorithms over triangular mesh algorithms such as accuracy, and computational resources have been extensively surveyed and discussed (Owen, 1998). So, it is both interesting and challenging to study and apply this type of algorithm to develop alternative FE methods to solve a 2-D magnetotelluric problem.
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
In this work, we applied the FE method with an alternative unstructured quadrilateral mesh to solve a 2-D magnetotelluric problem. An unstructured quadrilateral mesh was used for our FE codes and was generated using the paving algorithm (Blacker and Stephenson, 1991). Here, the paving algorithm is reviewed and summarised. Suitable FE codes with the selected mesh type were developed and tested. The resulting numerical solution was validated and compared to both an analytical solution and other numerical methods. Furthermore, the capabilities of the alternative mesh algorithm to handle the 2-D domain are presented and discussed. 2. Governing equations Magnetotelluric (MT) modeling describes the phenomena of the natural electromagnetic (EM) induction in the Earth. The naturally occurring EM fields act like a plane-wave with harmonic diffusion, the displacement currents can be neglected and the timedependence is assumed to be e −iyt , where y is the angular frequency. For the case where the electrical conductivity s is varied in only the x and z-directions, i.e., s = s(x, z), the second-order partial equations for the 2-D MT problem with the bounded region Y⊂R2 is given by
near complex boundaries of the domain and are unable to completely satisfy some constraints (Remacle et al., 2012). For these reasons, a direct approach was used for generating an unstructured quadrilateral mesh for the 2-D MT model employing the paving algorithm of Blacker and Stephenson (1991). An example of the paving algorithm for 2-D magnetotelluric model is illustrated in Fig. 1. At the beginning, the nodes are paved on the boundary of the domain (Fig. 1b). The distance between two nodes does not need to be uniform. Usually, the connectivities of nodes on the exterior C and interior Cint boundaries are set as permanent, i.e., the coordinate and number of node on these boundaries do not change. With permanent boundaries, an exterior paving boundary is paved counterclockwise and progresses inward from the exterior boundary. In contrast, the interior paving boundary is paved clockwise and progresses outward from the interior permanent boundary. This step is called row choice (Fig. 1c). Next, new rows are added and adjusted to correct the element size and improve the mesh quality by using
∂ ∂u ∂ ∂u a + a + bu = 0. ∂x ∂x ∂z ∂z
(1)
The notations a, b and u denote different representations depending on the two MT polarizations: E-polarization: u = Ey , H-polarization: u = Hy ,
a = 1,
b = iyls,
a = 1/s,
b = iyl,
(2) (3)
where Ey and Hy are the strike aligned electric and magnetic fields, respectively, l is the magnetic permeability in free space and l = l 0 = 4p × 10 −7 (V s/Am). The bounded region Y is defined by Y = Y1 ∪ Y2 ∪ C ∪ Cint , where Y1 is the air subregion, Y2 is the earth subregion, C is the outer boundary and Cint is the air-earth interface where electrical conductivity s is discontinuous. The partial differential Eq. (1) is subjected to the Dirichlet boundary conditions u = u0 (x, z)
on
C,
(4)
where u0 (x, z) is obtained by solving the one-dimensional MT problem. 3. Mesh algorithm To solve Eq. (1) for the 2-D magnetotelluric forward modeling subjected to the boundary condition as in Eq. (4) using the finite element method, discretization or mesh generation is the first and important step. For the unstructured quadrilateral mesh, a continuous 2-D domain, Y, was meshed into many subdomains or elements. e e Here, Y = ∪M e=1 Y , where Y is the e-th element with a quadrilateral shape and M is the total number of elements. Presently, there are two categories of generating unstructured quadrilateral mesh: direct and indirect approaches (Owen, 1998). The direct approach generates the mesh from the given domain whereas the indirect approach generates a quadrilateral mesh from an existing triangular mesh by splitting or merging elements. The algorithms of some indirect approaches provide elements of poor quality and do not guarantee complete transformation of triangular elements to quadrilateral elements (Owen, 1998). In contrast, direct approaches guarantee that the generated meshes are composed entirely of quadrilateral elements. However, some of them may produce low-quality elements
17
Fig. 1. Example of paving algorithm for 2-D model.
18
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
the row adjustment and smoothing processes, respectively (Fig. 1d). As rows begin to overlap or intersect at the interior of the subregion, they are connected together to form a valid quadrilateral mesh by using the seaming process (Fig. 1e). When the region has been completely paved, some local mesh is adjusted again to improve the local aspect ratios, reduce the number of irregular nodes, or eliminate elements with bad interior angles by deleting and/or adding
an element (Fig. 1f). Finally, the unstructured quadrilateral mesh is obtained. In this work, an unstructured quadrilateral mesh was generated using a paving algorithm of Trelis commercial software. The generated coordinates of all nodes and element-node data of all quadrilateral elements were required for our finite element codes developed in the course of the current study.
Fig. 2. The flowchart of finite element software.
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
19
4. Finite element codes In this section, the development of computer software for the FE method will be explained. The algorithm used to develop our software is illustrated in Fig. 2. For our software, the user must input the following necessary data: periods of wave, location of sites, structure of model defined by polygons, structure of boundary, information of all quadrilateral elements and nodes generated by paving algorithm. Then, the routine for reordering the data of elements and nodes using the Reverse Cuthill-McKee (RCM) algorithm is carried out to minimize the bandwidth of the resulting coefficient matrix and the nodes are classified as interior or exterior to be able to apply boundary conditions. Furthermore, the necessary information for solving the 1-D problem is set up and the resistivity is assigned to each element according to the model. The reordering process is required in this step because the efficiency of solving the system of equation with the direct method relies on the minimum of the bandwidth of the matrix. During the main process, the coefficient matrix for each element is created and assembled to the augmented coefficient matrix by assembling all elemental matrices and applying the Dirichlet boundary conditions obtained by solving the 1-D problem, the system of equations is finally expressed as Ku = b,
(5)
where the coefficient matrix K is large, sparse, symmetric and complex but not Hermitian. u is the unknown vector and b is the source vector corresponding to the boundary conditions. The system of equations is solved by LU factorization using partial pivoting with row interchange. In the final process, the interpolation of the orthogonal field and the responses including impedance, apparent resistivity and phase at each site are estimated. These responses are the output of our software. Note that all subroutines described here were developed in MATLAB without MATLAB’s special toolboxes. The codes can be run on Windows, OS X or Linux systems with a MATLAB installation. 5. Numerical results In this section, the efficiency, accuracy, and reliability of our alternative FE codes are presented, compared and discussed. Firstly, the half-space model was used as the test model to present its accuracy and this model was also used for comparison with another FE method. The Root Mean Square (RMS) error was used as a measure of the accuracy of the calculated MT responses. For the reliability of our codes, the COMMEMI2D-4 model presented by Zhdanov et al. (1997) is used as the test model. Finally, the model of the sea is used to present its flexibility. Furthermore, simple ocean coast effects upon magnetotelluric responses are presented and discussed. 5.1. Testing accuracy The test of the accuracy of our FE method began with a half-space model. The half-space model is illustrated in Fig. 3 (top). The size of MT domain was designed as 200 km × 200 km and the thickness of the air layer extends to 200 km. The earth’s surface was assumed to be flat and homogeneous with a resistivity of 100 ohm-m. We supposed that there were 19 MT sites with a distance of 3 km located on the earth’s surface (z = 0 km) in the range −30 km to 30 km. The MT responses including the apparent resistivity and phases were computed at 4 periods: T = 0.1, 1.0, 10, and 100 s for both E- and H-polarizations. To solve the 2-D MT problem using the finite element method with an unstructured quadrilateral mesh, the given domain was discretized into many elements and nodes using a paving algorithm
Fig. 3. The 100ohm-m homogeneous earth (top) and meshed with unstructured quadrilateral elements (bottom). The red triangles represent the locations of MT sites.
from the Trelis mesh software. The computing domain and its generated mesh are shown in Fig. 3 (bottom). The smallest element size (minimum edge length of an element) around the MT sites was set to 150 m. The size of the elements located far from the sites and the area of interest were increased at a rate of 1.2 and 1.4 for the earth and air layers, respectively. The total numbers of nodes and elements generated were 16,885 and 16,848, respectively. Table 1 RMS errors in apparent resistivity and phase. Period (s)
0.1 1.0 10 100
E-polarization
H-polarization
App. res.
Phase
App. res.
Phase
3.9194 1.3615 0.4068 0.1659
1.0714 0.3908 0.1238 0.0583
3.6532 1.3284 0.4066 0.14316
1.0644 0.3878 0.1240 0.0574
20
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24 Table 2 Mesh data of two numerical methods: FET and FEQ methods. Mesh type
Coarse Medium Fine
Smallest element
Biggest element
500 250 100
10,000 10,000 10,000
Increasing rate
1.23 1.1 1.1
Num. of elements
Num. of nodes
FET
FEQ
FET
FEQ
7528 15,270 32,960
3383 7207 16,616
3810 7697 16,567
3440 7264 16,729
Table 3 RMS error of MT responses at T = 1s. Mesh types
E-polarization
H-polarization
FET
Coarse Medium Fine
FEQ
FET
FEQ
qxy
0xy
qyx
0yx
qyx
0yx
qyx
0yx
14.7279 10.0093 9.4070
3.7370 2.6505 2.5047
12.5100 6.7424 5.3256
3.1222 1.7817 1.4288
12.9469 9.2216 8.4740
3.7692 2.6938 2.4884
10.5608 6.1388 4.9774
3.1256 1.7808 1.4347
For each mesh type, the numbers of nodes used for the two FE methods were set as close as possible since these parameters affect the accuracy of MT responses. To compare the accuracy of the two methods, the RMS errors of MT responses were calculated using two methods. The 31 MT sites were set in the range −30 km to 30 km with a distance of 2 km and their MT responses were computed at periods T = 1, 3.1623, and 10 s, respectively. The RMS errors at T = 1s, are shown in Table 3. The maximum RMS errors of the two methods occurred using a coarse mesh whereas the minimum occurred using a fine mesh, i.e., the RMS errors were reduced when the mesh size was reduced. All RMS errors for FET were greater than those of FEQ for all mesh types. This indicates that the FEQ method was more accurate than the FET method. RMS errors at T = 3.1623s, are shown in Table 4. As in the previous case, the maximum RMS errors of the two methods occurred when using a coarse mesh and the smallest were observed with the fine mesh. However, the levels of all RMS errors were lower than those of the previous period. All RMS errors generated by FET were greater than those of FEQ for all mesh types. For the last period, the RMS errors of the last period, T = 10s, are shown in Table 5. As with the previous periods, the maximum RMS errors of the two methods occurred with the coarse mesh. However, the levels of all RMS errors were the lowest when compared to those of all previous periods. The RMS errors were reduced when the mesh
For each period, our FE algorithm took about 34 s to complete the 2-D MT forward modeling computations. For this model, the exact apparent resistivity for each MT site was identical to the real resistivity, i.e., 100 ohm-m and the phase was equal to 45◦ . The RMS errors of the apparent resistivity and phase are shown in Table 1. It is obvious that the RMS error was reduced when the period was increased, with the maximum error at the shortest period, T = 0.1s, and the minimum at the longest period, T = 100s. The results indicate that our FE codes provide very good accuracy. 5.2. Comparison of accuracy and efficiency In this section, the efficiency and accuracy of our FE algorithm (FEQ) are compared to another FE method with a structured triangular mesh (FET). To show a comparison of accuracy, the test was again based on a 100 ohm-m half-space model. The model was meshed at three levels: coarse, medium and fine mesh. The meshing algorithm for each numerical method was controlled using three parameters including the smallest element size (in meter), rate of increase of element size and the biggest element size (in meter). Using those parameters, the number of nodes and elements were optimized by the mesh algorithm. It is notable that, the unstructured triangular meshes were generated using the Delaunay algorithm from the Trelis mesh software. The characteristics of the generated meshes are summarised in Table 2. Table 4 RMS error of MT responses at T = 3.1623s. Mesh types
E-polarization
H-polarization
FET
Coarse Medium Fine
FEQ
FET
FEQ
qxy
0xy
qyx
0yx
qyx
0yx
qyx
0yx
7.6519 5.0858 4.7836
2.0630 1.4536 1.3949
6.8111 3.7381 2.9590
1.7930 1.0134 0.8110
8.0460 5.6902 5.1700
2.1115 1.4605 1.3487
6.1758 3.5369 2.8560
1.7938 1.0136 0.8166
Table 5 RMS error of MT responses at T = 10s. Mesh types
E-polarization
H-polarization
FET
Coarse Medium Fine
FEQ
FET
FEQ
qxy
0xy
qyx
0yx
qyx
0yx
qyx
0yx
4.8468 3.4762 3.1693
1.0869 0.9020 0.9251
3.7754 2.1023 1.6627
1.0228 0.5743 0.4593
6.4328 5.4320 5.1236
2.1540 1.6909 1.5764
3.5582 2.0284 1.6285
1.0209 0.5723 0.4619
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
Fig. 4. Comparison of the usage of CPU time of FET (blue) and FEQ (green) methods versus different mesh types.
Fig. 5. The COMMEMI2D-4 model with the location of MT sites (red triangle) and its meshing near air-earth’s interface and at shallow depths.
21
size was reduced and all RMS errors calculated by FET were still greater than those of FEQ for all mesh types. In these experiments, our FEQ method showed better accuracy than the FET method. The accuracy of FET method can be improved using an adaptive scheme (Key and Weiss, 2006; Franke et al., 2007). However, the adaptive scheme generates a huge number of nodes and elements during mesh generation thereby generating higher computational overhead and it may take a longer calculation time. Next, the efficiency of our FEQ method in term of total CPU time was compared to FET. The usage of CPU time per period for the two numerical methods is plotted versus mesh type and shown in Fig. 4. The total CPU times of the two methods increased when the mesh type was changed from coarse to fine mesh because the size of the system of equations was increased. The FEQ method construct an elemental matrix with dimensions of 4 × 4 and its members were evaluated by using 4-point Gauss-quadrature numerical integration, whereas the FET constructs an elemental matrix with dimensions of 3 × 3 and its members were evaluated directly without using any
Fig. 6. Apparent resistivity at period T = 1 s for E-polarization (top) and Hpolarization (bottom).
22
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
numerical integration. The augmented coefficient matrix for all elements in the matrix of the FEQ method is sparse as is that of the FET method, but its nonzero entries were denser. For these reasons, our FEQ method was slower than the FET method. 5.3. Reliability In this section, the reliability of our FEQ codes is presented and confirmed by testing on the COMMEMI2D-4 model. The COMMEMI2D-4 model, proposed by Zhdanov et al. (1997) for comparing of modeling methods for electromagnetic induction, is illustrated in Fig. 5 (top). The resistivity structure of the model consists of three layers with resistivity 25 ohm-m at the top, 1000 ohm-m in the middle and 5 ohm-m at the bottom of the model. The embedded 2.5 ohm-m graben-like structure tapers off as a layer to the right-hand of the model boundary. On the left-hand model boundary, the stratification was different because of a 10 ohm-m semi-infinite
Fig. 7. Apparent resistivity at period T = 100 s for E-polarization (top) and Hpolarization (bottom).
layer. For this test, the size of the model was set to 200 km × 200 km. To do the 2-D MT modeling, the model was discretized into 18,990 nodes and 18,910 elements. The smallest element size was 80 m. Some mesh points near earth’s surface and at shallow depths with local refinement of the mesh are shown in Fig. 5 (bottom). The apparent resistivity of 23 sites, located in the range −30 km to 30 km., were computed at periods T = 1s and 100 s. Fig. 6 shows the apparent resistivity values at period T = 1s, while Fig. 7 gives those values at T = 100s. Our results were quite comparable with those of the COMMEMI2D-4 projects. The errors of the apparent resistivity of the E-polarization were less than those of the H-polarization. These comparisons confirm the validity and reliability of our FEQ codes. 5.4. Flexibility The flexibility of our codes and the example of ocean coast effects on magnetotelluric responses are presented and discussed in this section. Initially, we designed the model as a composition of the
Fig. 8. The coast model (top) and its mesh (bottom). The red triangles represent the location of MT sites.
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
ocean and earth as shown in Fig. 8 (top). The ocean’s level and earth’s surface are located at z = 0 km. The coast is located at x = −10 km. The seafloor was varied from shallow to deep. The deepest of part of the ocean is at the left boundary with a depth of 5 km. For the computation, the resistivity of the ocean and earth were set as 0.2 ohm-m and 100 ohm-m, respectively. The 21 onshore sites were located in the range of −9 km to 50 km with a distance of 3 km and the apparent resistivity and phase were calculated at periods of T = 0.1, 1, 10, 100 s for each site. The mesh of the model included 16,839 nodes and 16,788 elements as shown in Fig. 8 (bottom). Local refinement of the mesh was applied near the earth’s surface where MT sites were located to improve the accuracy of the responses. The obtained apparent resistivity and phase are shown in Figs. 9 and 10, respectively. The results showed that the largest effect appeared at MT sites nearest to the ocean, whereas the effect was the smallest at the end site. Additionally, the coast effect dominated the responses at longer periods, whereas it had no effect at shorter periods. For E-polarization, the apparent resistivity at sites located near the coast was reduced compared to those that were distant. In contrast, the phase increased at nearby sites. For H-polarization,
Fig. 9. The apparent resistivity (log scale) for E-polarization (top) and H-polarization (bottom).
23
the coast effect of the apparent resistivity was the same as for Epolarization, but it was more dominant. Additionally, the effect of the coast on the phase was opposite to that of E-polarization. This model is an example of a study of the coast effects using magnetotelluric data that relies on 2-D magnetotelluric modeling employing our alternative FE method. Previous studies of the coast effects of magnetotelluric data have been presented in the literature (Malleswari and Veeraswamy, 2014; Pandey et al., 2008; Singh et al., 1995). 6. Discussions and conclusions The finite element method with an unstructured quadrilateral mesh was successfully developed and applied to solve a 2-D magnetotelluric problem. Our FE codes gave good accuracy of response when it was tested on a simple half-space model. With the appropriate mesh, the maximum RMS of obtained responses was below 4% for all test periods. Additionally, our codes gave better accuracy of response when compared to those of the FE method with unstructured triangular meshes for all test periods and mesh types. When the mesh type was changed from a coarser to a finer mesh, the maximum RMS was reduced from about 14% to 5% at the shortest period and from about 4% to 2% at the longest period. However, our FE codes
Fig. 10. The phase for E-polarization (top) and H-polarization (bottom).
24
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
in terms of CPU time are slower than those of the FE method with an unstructured triangular mesh because the processes to evaluate the elemental matrix and coefficient matrix are more complicated. When evaluating the reliability of the codes, it was found that our responses and those calculated by the COMMEMI projects (Zhdanov et al., 1997) are comparable for both polarizations. This shows that our codes are accurate and reliable. Furthermore, our FE codes using an unstructured quadrilateral mesh provide additional power features including local mesh refinement and the capability of handling complex structures, which are useful for a complex 2-D domain. This indicates that our FE codes are flexible and appropriate for use in any real 2-D domain. In the future, these forward modeling codes can be used as a forward calculator to further develop 2-D magnetotelluric inversion for interpreting real magnetotelluric data. Acknowledgments This work was supported by the DPST Research grant [grant number 002/2557].
References Bern, M., 1997. Quadrilateral meshing by circle packing. Int. J. Comput. Geom. Appl. 7–20. Blacker, T.D., Stephenson, M.B., 1991. Paving. A new approach to automated quadrilateral mesh generation. Int. J. Numer. Methods Eng. 32 (4), 811–847.
Franke, A., Börner, R.U., Spitzer, K., 2007. Adaptive unstructured grid finite element simulation of two-dimensional magnetotelluric fields for arbitrary surface and seafloor topography. Geophys. J. Int. 171 (1), 71–86. Key, K., Weiss, C., 2006. Adaptive finite-element modeling using unstructured grids: the 2d magnetotelluric example. Geophysics 71 (6), G291–G299. Lee, S.K., Kim, H., Song, Y., Lee, C.K., 2009. MT2DInvmatlab-a program in MATLAB and FORTRAN for two-dimensional magnetotelluric inversion. Comput. Geosci. 35 (8), 1722–1734. Malleswari, D., Veeraswamy, K., 2014. Numerical simulation of coast effect on magnetotelluric measurements. Acta Geodaet. Geophys. 49 (1), 17–35. Owen, S.J., 1998. A Survey of Unstructured Mesh Generation Technology. International Meshing Roundtable pp. 239–267 Pandey, D., Sinha, M., MacGregor, L., Singh, S., 2008. Ocean coast effect on magnetotelluric data: a case study from Kachchh, India. Mar. Geophys. Res. 29 (3), 185–193. Rao, K.P., Babu, G.A., 2006. EMOD2d–a program in C++ for finite difference modelling of magnetotelluric TM mode responses over 2D earth. Comput. Geosci. 32 (9), 1499–1511. Reddy, I.K., Rankin, D., 1975. Magnetotelluric response of laterally inhomogeneous and anisotropic media. Geophysics 40 (6), 1035–1045. Remacle, J.F., Lambrechts, J., Seny, B., Marchandise, E., Johnen, A., Geuzainet, C., 2012. Blossom-quad: a non-uniform quadrilateral mesh generator using a minimum-cost perfect-matching algorithm. Int. J. Numer. Methods Eng. 89 (9), 1102–1119. Sarrate, J., Huerta, A., 2000. Efficient unstructured quadrilateral mesh generation. Int. J. Numer. Methods Eng. 49 (10), 1327–1350. Singh, U.K., Kant, Y., Singh, R.P., 1995. Effect of coast on magnetotelluric measurements in India. Ann. Geophys. 38 (3-4). Wannamaker, P.E., Stodt, J.A., Rijo, L., 1986. Two-dimensional topographic responses in magnetotellurics modeled using finite elements. Geophysics 51 (11), 2131–2144. Zhdanov, M.S., Varentsov, I.M., Weaver, J.T., Golubev, N.G., Krylov, V.A., 1997. Methods for modelling electromagnetic fields results from COMMEMI - the international project on the comparison of modelling methods for electromagnetic induction. J. Appl. Geophys. 37 (3-4), 133–271.
Contents lists available at ScienceDirect
Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo
2-D magnetotelluric modeling using finite element method incorporating unstructured quadrilateral elements Weerachai Sarakorn Department of Mathematics, Faculty of Science, Khon Kaen University, Khon Kaen, Thailand
A R T I C L E
I N F O
Article history: Received 25 May 2016 Received in revised form 10 February 2017 Accepted 13 February 2017 Available online 16 February 2017 JEL classification: 02.70.Dh 91.25. − r 2000 MSC: 65N30 65N50 74S05
A B S T R A C T In this research, the finite element (FE) method incorporating quadrilateral elements for solving 2-D MT modeling was presented. The finite element software was developed, employing a paving algorithm to generate the unstructured quadrilateral mesh. The accuracy, efficiency, reliability, and flexibility of our FE forward modeling are presented, compared and discussed. The numerical results indicate that our FE codes using an unstructured quadrilateral mesh provide good accuracy when the local mesh refinement is applied around sites and in the area of interest, with superior results when compared to other FE methods. The reliability of the developed codes was also confirmed when comparing both analytical solutions and COMMEMI2D model. Furthermore, our developed FE codes incorporating an unstructured quadrilateral mesh showed useful and powerful features such as handling irregular and complex subregions and providing local refinement of the mesh for a 2-D domain as closely as unstructured triangular mesh but it requires less number of elements in a mesh. © 2017 Elsevier B.V. All rights reserved.
Keywords: Magnetotelluric modeling Finite element method Unstructured quadrilateral mesh Paving algorithm
1. Introduction From the past to the present, two-dimensional magnetotelluric modeling was successfully solved by various numerical techniques such as the finite difference (FD) method (Rao and Babu, 2006) and finite element (FE) method (Reddy and Rankin, 1975; Wannamaker et al., 1986; Key and Weiss, 2006; Franke et al., 2007; Lee et al., 2009). The efficiency and accuracy of the FD method were proven for a simple 2-D model (Rao and Babu, 2006). With the limitation of rectangular grids, many complex subregions in the 2-D domain such as modeling topography, irregular anomalous bodies, and seafloor cannot be handled by its discretization. Therefore, the solution obtained by the FD method may be inaccurate and unrelialable. The finite element (FE) method is another numerical technique that is often used for 2-D MT forward modeling. Its efficiency and accuracy have been proven (Reddy and Rankin, 1975; Wannamaker et al., 1986; Key and Weiss, 2006; Franke et al., 2007; Lee et al., 2009). Unlike the FD method, the FE method is more flexible on discretization. Various mesh schemes such as structured triangular meshes (Wannamaker et al., 1986), unstructured triangular meshes with
http://dx.doi.org/10.1016/j.jappgeo.2017.02.005 0009-2541/© 2017 Elsevier B.V. All rights reserved.
adaptive refinement (Franke et al., 2007; Key and Weiss, 2006) and structured quadrilateral meshes (Reddy and Rankin, 1975; Lee et al., 2009) can be used in the FE method. The structured triangular and quadrilateral mesh algorithms are simple and can handle domains with complex subregions, but the local refinement mesh feature cannot be implemented efficiently. In contrast, the unstructured triangular mesh algorithm can overcome this problem. However, it needs to incorporate an adaptive scheme to improve its accuracy and it requires a large number of triangular elements in complex subregions (Key and Weiss, 2006; Franke et al., 2007). Presently, implementations of unstructured quadrilateral meshes already in the FE method are popular as discussed by Blacker and Stephenson (1991); Bern (1997); Owen (1998); Sarrate and Huerta (2000); Remacle et al. (2012). The advantages of unstructured quadrilateral mesh algorithms over triangular mesh algorithms such as accuracy, and computational resources have been extensively surveyed and discussed (Owen, 1998). So, it is both interesting and challenging to study and apply this type of algorithm to develop alternative FE methods to solve a 2-D magnetotelluric problem.
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
In this work, we applied the FE method with an alternative unstructured quadrilateral mesh to solve a 2-D magnetotelluric problem. An unstructured quadrilateral mesh was used for our FE codes and was generated using the paving algorithm (Blacker and Stephenson, 1991). Here, the paving algorithm is reviewed and summarised. Suitable FE codes with the selected mesh type were developed and tested. The resulting numerical solution was validated and compared to both an analytical solution and other numerical methods. Furthermore, the capabilities of the alternative mesh algorithm to handle the 2-D domain are presented and discussed. 2. Governing equations Magnetotelluric (MT) modeling describes the phenomena of the natural electromagnetic (EM) induction in the Earth. The naturally occurring EM fields act like a plane-wave with harmonic diffusion, the displacement currents can be neglected and the timedependence is assumed to be e −iyt , where y is the angular frequency. For the case where the electrical conductivity s is varied in only the x and z-directions, i.e., s = s(x, z), the second-order partial equations for the 2-D MT problem with the bounded region Y⊂R2 is given by
near complex boundaries of the domain and are unable to completely satisfy some constraints (Remacle et al., 2012). For these reasons, a direct approach was used for generating an unstructured quadrilateral mesh for the 2-D MT model employing the paving algorithm of Blacker and Stephenson (1991). An example of the paving algorithm for 2-D magnetotelluric model is illustrated in Fig. 1. At the beginning, the nodes are paved on the boundary of the domain (Fig. 1b). The distance between two nodes does not need to be uniform. Usually, the connectivities of nodes on the exterior C and interior Cint boundaries are set as permanent, i.e., the coordinate and number of node on these boundaries do not change. With permanent boundaries, an exterior paving boundary is paved counterclockwise and progresses inward from the exterior boundary. In contrast, the interior paving boundary is paved clockwise and progresses outward from the interior permanent boundary. This step is called row choice (Fig. 1c). Next, new rows are added and adjusted to correct the element size and improve the mesh quality by using
∂ ∂u ∂ ∂u a + a + bu = 0. ∂x ∂x ∂z ∂z
(1)
The notations a, b and u denote different representations depending on the two MT polarizations: E-polarization: u = Ey , H-polarization: u = Hy ,
a = 1,
b = iyls,
a = 1/s,
b = iyl,
(2) (3)
where Ey and Hy are the strike aligned electric and magnetic fields, respectively, l is the magnetic permeability in free space and l = l 0 = 4p × 10 −7 (V s/Am). The bounded region Y is defined by Y = Y1 ∪ Y2 ∪ C ∪ Cint , where Y1 is the air subregion, Y2 is the earth subregion, C is the outer boundary and Cint is the air-earth interface where electrical conductivity s is discontinuous. The partial differential Eq. (1) is subjected to the Dirichlet boundary conditions u = u0 (x, z)
on
C,
(4)
where u0 (x, z) is obtained by solving the one-dimensional MT problem. 3. Mesh algorithm To solve Eq. (1) for the 2-D magnetotelluric forward modeling subjected to the boundary condition as in Eq. (4) using the finite element method, discretization or mesh generation is the first and important step. For the unstructured quadrilateral mesh, a continuous 2-D domain, Y, was meshed into many subdomains or elements. e e Here, Y = ∪M e=1 Y , where Y is the e-th element with a quadrilateral shape and M is the total number of elements. Presently, there are two categories of generating unstructured quadrilateral mesh: direct and indirect approaches (Owen, 1998). The direct approach generates the mesh from the given domain whereas the indirect approach generates a quadrilateral mesh from an existing triangular mesh by splitting or merging elements. The algorithms of some indirect approaches provide elements of poor quality and do not guarantee complete transformation of triangular elements to quadrilateral elements (Owen, 1998). In contrast, direct approaches guarantee that the generated meshes are composed entirely of quadrilateral elements. However, some of them may produce low-quality elements
17
Fig. 1. Example of paving algorithm for 2-D model.
18
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
the row adjustment and smoothing processes, respectively (Fig. 1d). As rows begin to overlap or intersect at the interior of the subregion, they are connected together to form a valid quadrilateral mesh by using the seaming process (Fig. 1e). When the region has been completely paved, some local mesh is adjusted again to improve the local aspect ratios, reduce the number of irregular nodes, or eliminate elements with bad interior angles by deleting and/or adding
an element (Fig. 1f). Finally, the unstructured quadrilateral mesh is obtained. In this work, an unstructured quadrilateral mesh was generated using a paving algorithm of Trelis commercial software. The generated coordinates of all nodes and element-node data of all quadrilateral elements were required for our finite element codes developed in the course of the current study.
Fig. 2. The flowchart of finite element software.
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
19
4. Finite element codes In this section, the development of computer software for the FE method will be explained. The algorithm used to develop our software is illustrated in Fig. 2. For our software, the user must input the following necessary data: periods of wave, location of sites, structure of model defined by polygons, structure of boundary, information of all quadrilateral elements and nodes generated by paving algorithm. Then, the routine for reordering the data of elements and nodes using the Reverse Cuthill-McKee (RCM) algorithm is carried out to minimize the bandwidth of the resulting coefficient matrix and the nodes are classified as interior or exterior to be able to apply boundary conditions. Furthermore, the necessary information for solving the 1-D problem is set up and the resistivity is assigned to each element according to the model. The reordering process is required in this step because the efficiency of solving the system of equation with the direct method relies on the minimum of the bandwidth of the matrix. During the main process, the coefficient matrix for each element is created and assembled to the augmented coefficient matrix by assembling all elemental matrices and applying the Dirichlet boundary conditions obtained by solving the 1-D problem, the system of equations is finally expressed as Ku = b,
(5)
where the coefficient matrix K is large, sparse, symmetric and complex but not Hermitian. u is the unknown vector and b is the source vector corresponding to the boundary conditions. The system of equations is solved by LU factorization using partial pivoting with row interchange. In the final process, the interpolation of the orthogonal field and the responses including impedance, apparent resistivity and phase at each site are estimated. These responses are the output of our software. Note that all subroutines described here were developed in MATLAB without MATLAB’s special toolboxes. The codes can be run on Windows, OS X or Linux systems with a MATLAB installation. 5. Numerical results In this section, the efficiency, accuracy, and reliability of our alternative FE codes are presented, compared and discussed. Firstly, the half-space model was used as the test model to present its accuracy and this model was also used for comparison with another FE method. The Root Mean Square (RMS) error was used as a measure of the accuracy of the calculated MT responses. For the reliability of our codes, the COMMEMI2D-4 model presented by Zhdanov et al. (1997) is used as the test model. Finally, the model of the sea is used to present its flexibility. Furthermore, simple ocean coast effects upon magnetotelluric responses are presented and discussed. 5.1. Testing accuracy The test of the accuracy of our FE method began with a half-space model. The half-space model is illustrated in Fig. 3 (top). The size of MT domain was designed as 200 km × 200 km and the thickness of the air layer extends to 200 km. The earth’s surface was assumed to be flat and homogeneous with a resistivity of 100 ohm-m. We supposed that there were 19 MT sites with a distance of 3 km located on the earth’s surface (z = 0 km) in the range −30 km to 30 km. The MT responses including the apparent resistivity and phases were computed at 4 periods: T = 0.1, 1.0, 10, and 100 s for both E- and H-polarizations. To solve the 2-D MT problem using the finite element method with an unstructured quadrilateral mesh, the given domain was discretized into many elements and nodes using a paving algorithm
Fig. 3. The 100ohm-m homogeneous earth (top) and meshed with unstructured quadrilateral elements (bottom). The red triangles represent the locations of MT sites.
from the Trelis mesh software. The computing domain and its generated mesh are shown in Fig. 3 (bottom). The smallest element size (minimum edge length of an element) around the MT sites was set to 150 m. The size of the elements located far from the sites and the area of interest were increased at a rate of 1.2 and 1.4 for the earth and air layers, respectively. The total numbers of nodes and elements generated were 16,885 and 16,848, respectively. Table 1 RMS errors in apparent resistivity and phase. Period (s)
0.1 1.0 10 100
E-polarization
H-polarization
App. res.
Phase
App. res.
Phase
3.9194 1.3615 0.4068 0.1659
1.0714 0.3908 0.1238 0.0583
3.6532 1.3284 0.4066 0.14316
1.0644 0.3878 0.1240 0.0574
20
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24 Table 2 Mesh data of two numerical methods: FET and FEQ methods. Mesh type
Coarse Medium Fine
Smallest element
Biggest element
500 250 100
10,000 10,000 10,000
Increasing rate
1.23 1.1 1.1
Num. of elements
Num. of nodes
FET
FEQ
FET
FEQ
7528 15,270 32,960
3383 7207 16,616
3810 7697 16,567
3440 7264 16,729
Table 3 RMS error of MT responses at T = 1s. Mesh types
E-polarization
H-polarization
FET
Coarse Medium Fine
FEQ
FET
FEQ
qxy
0xy
qyx
0yx
qyx
0yx
qyx
0yx
14.7279 10.0093 9.4070
3.7370 2.6505 2.5047
12.5100 6.7424 5.3256
3.1222 1.7817 1.4288
12.9469 9.2216 8.4740
3.7692 2.6938 2.4884
10.5608 6.1388 4.9774
3.1256 1.7808 1.4347
For each mesh type, the numbers of nodes used for the two FE methods were set as close as possible since these parameters affect the accuracy of MT responses. To compare the accuracy of the two methods, the RMS errors of MT responses were calculated using two methods. The 31 MT sites were set in the range −30 km to 30 km with a distance of 2 km and their MT responses were computed at periods T = 1, 3.1623, and 10 s, respectively. The RMS errors at T = 1s, are shown in Table 3. The maximum RMS errors of the two methods occurred using a coarse mesh whereas the minimum occurred using a fine mesh, i.e., the RMS errors were reduced when the mesh size was reduced. All RMS errors for FET were greater than those of FEQ for all mesh types. This indicates that the FEQ method was more accurate than the FET method. RMS errors at T = 3.1623s, are shown in Table 4. As in the previous case, the maximum RMS errors of the two methods occurred when using a coarse mesh and the smallest were observed with the fine mesh. However, the levels of all RMS errors were lower than those of the previous period. All RMS errors generated by FET were greater than those of FEQ for all mesh types. For the last period, the RMS errors of the last period, T = 10s, are shown in Table 5. As with the previous periods, the maximum RMS errors of the two methods occurred with the coarse mesh. However, the levels of all RMS errors were the lowest when compared to those of all previous periods. The RMS errors were reduced when the mesh
For each period, our FE algorithm took about 34 s to complete the 2-D MT forward modeling computations. For this model, the exact apparent resistivity for each MT site was identical to the real resistivity, i.e., 100 ohm-m and the phase was equal to 45◦ . The RMS errors of the apparent resistivity and phase are shown in Table 1. It is obvious that the RMS error was reduced when the period was increased, with the maximum error at the shortest period, T = 0.1s, and the minimum at the longest period, T = 100s. The results indicate that our FE codes provide very good accuracy. 5.2. Comparison of accuracy and efficiency In this section, the efficiency and accuracy of our FE algorithm (FEQ) are compared to another FE method with a structured triangular mesh (FET). To show a comparison of accuracy, the test was again based on a 100 ohm-m half-space model. The model was meshed at three levels: coarse, medium and fine mesh. The meshing algorithm for each numerical method was controlled using three parameters including the smallest element size (in meter), rate of increase of element size and the biggest element size (in meter). Using those parameters, the number of nodes and elements were optimized by the mesh algorithm. It is notable that, the unstructured triangular meshes were generated using the Delaunay algorithm from the Trelis mesh software. The characteristics of the generated meshes are summarised in Table 2. Table 4 RMS error of MT responses at T = 3.1623s. Mesh types
E-polarization
H-polarization
FET
Coarse Medium Fine
FEQ
FET
FEQ
qxy
0xy
qyx
0yx
qyx
0yx
qyx
0yx
7.6519 5.0858 4.7836
2.0630 1.4536 1.3949
6.8111 3.7381 2.9590
1.7930 1.0134 0.8110
8.0460 5.6902 5.1700
2.1115 1.4605 1.3487
6.1758 3.5369 2.8560
1.7938 1.0136 0.8166
Table 5 RMS error of MT responses at T = 10s. Mesh types
E-polarization
H-polarization
FET
Coarse Medium Fine
FEQ
FET
FEQ
qxy
0xy
qyx
0yx
qyx
0yx
qyx
0yx
4.8468 3.4762 3.1693
1.0869 0.9020 0.9251
3.7754 2.1023 1.6627
1.0228 0.5743 0.4593
6.4328 5.4320 5.1236
2.1540 1.6909 1.5764
3.5582 2.0284 1.6285
1.0209 0.5723 0.4619
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
Fig. 4. Comparison of the usage of CPU time of FET (blue) and FEQ (green) methods versus different mesh types.
Fig. 5. The COMMEMI2D-4 model with the location of MT sites (red triangle) and its meshing near air-earth’s interface and at shallow depths.
21
size was reduced and all RMS errors calculated by FET were still greater than those of FEQ for all mesh types. In these experiments, our FEQ method showed better accuracy than the FET method. The accuracy of FET method can be improved using an adaptive scheme (Key and Weiss, 2006; Franke et al., 2007). However, the adaptive scheme generates a huge number of nodes and elements during mesh generation thereby generating higher computational overhead and it may take a longer calculation time. Next, the efficiency of our FEQ method in term of total CPU time was compared to FET. The usage of CPU time per period for the two numerical methods is plotted versus mesh type and shown in Fig. 4. The total CPU times of the two methods increased when the mesh type was changed from coarse to fine mesh because the size of the system of equations was increased. The FEQ method construct an elemental matrix with dimensions of 4 × 4 and its members were evaluated by using 4-point Gauss-quadrature numerical integration, whereas the FET constructs an elemental matrix with dimensions of 3 × 3 and its members were evaluated directly without using any
Fig. 6. Apparent resistivity at period T = 1 s for E-polarization (top) and Hpolarization (bottom).
22
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
numerical integration. The augmented coefficient matrix for all elements in the matrix of the FEQ method is sparse as is that of the FET method, but its nonzero entries were denser. For these reasons, our FEQ method was slower than the FET method. 5.3. Reliability In this section, the reliability of our FEQ codes is presented and confirmed by testing on the COMMEMI2D-4 model. The COMMEMI2D-4 model, proposed by Zhdanov et al. (1997) for comparing of modeling methods for electromagnetic induction, is illustrated in Fig. 5 (top). The resistivity structure of the model consists of three layers with resistivity 25 ohm-m at the top, 1000 ohm-m in the middle and 5 ohm-m at the bottom of the model. The embedded 2.5 ohm-m graben-like structure tapers off as a layer to the right-hand of the model boundary. On the left-hand model boundary, the stratification was different because of a 10 ohm-m semi-infinite
Fig. 7. Apparent resistivity at period T = 100 s for E-polarization (top) and Hpolarization (bottom).
layer. For this test, the size of the model was set to 200 km × 200 km. To do the 2-D MT modeling, the model was discretized into 18,990 nodes and 18,910 elements. The smallest element size was 80 m. Some mesh points near earth’s surface and at shallow depths with local refinement of the mesh are shown in Fig. 5 (bottom). The apparent resistivity of 23 sites, located in the range −30 km to 30 km., were computed at periods T = 1s and 100 s. Fig. 6 shows the apparent resistivity values at period T = 1s, while Fig. 7 gives those values at T = 100s. Our results were quite comparable with those of the COMMEMI2D-4 projects. The errors of the apparent resistivity of the E-polarization were less than those of the H-polarization. These comparisons confirm the validity and reliability of our FEQ codes. 5.4. Flexibility The flexibility of our codes and the example of ocean coast effects on magnetotelluric responses are presented and discussed in this section. Initially, we designed the model as a composition of the
Fig. 8. The coast model (top) and its mesh (bottom). The red triangles represent the location of MT sites.
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
ocean and earth as shown in Fig. 8 (top). The ocean’s level and earth’s surface are located at z = 0 km. The coast is located at x = −10 km. The seafloor was varied from shallow to deep. The deepest of part of the ocean is at the left boundary with a depth of 5 km. For the computation, the resistivity of the ocean and earth were set as 0.2 ohm-m and 100 ohm-m, respectively. The 21 onshore sites were located in the range of −9 km to 50 km with a distance of 3 km and the apparent resistivity and phase were calculated at periods of T = 0.1, 1, 10, 100 s for each site. The mesh of the model included 16,839 nodes and 16,788 elements as shown in Fig. 8 (bottom). Local refinement of the mesh was applied near the earth’s surface where MT sites were located to improve the accuracy of the responses. The obtained apparent resistivity and phase are shown in Figs. 9 and 10, respectively. The results showed that the largest effect appeared at MT sites nearest to the ocean, whereas the effect was the smallest at the end site. Additionally, the coast effect dominated the responses at longer periods, whereas it had no effect at shorter periods. For E-polarization, the apparent resistivity at sites located near the coast was reduced compared to those that were distant. In contrast, the phase increased at nearby sites. For H-polarization,
Fig. 9. The apparent resistivity (log scale) for E-polarization (top) and H-polarization (bottom).
23
the coast effect of the apparent resistivity was the same as for Epolarization, but it was more dominant. Additionally, the effect of the coast on the phase was opposite to that of E-polarization. This model is an example of a study of the coast effects using magnetotelluric data that relies on 2-D magnetotelluric modeling employing our alternative FE method. Previous studies of the coast effects of magnetotelluric data have been presented in the literature (Malleswari and Veeraswamy, 2014; Pandey et al., 2008; Singh et al., 1995). 6. Discussions and conclusions The finite element method with an unstructured quadrilateral mesh was successfully developed and applied to solve a 2-D magnetotelluric problem. Our FE codes gave good accuracy of response when it was tested on a simple half-space model. With the appropriate mesh, the maximum RMS of obtained responses was below 4% for all test periods. Additionally, our codes gave better accuracy of response when compared to those of the FE method with unstructured triangular meshes for all test periods and mesh types. When the mesh type was changed from a coarser to a finer mesh, the maximum RMS was reduced from about 14% to 5% at the shortest period and from about 4% to 2% at the longest period. However, our FE codes
Fig. 10. The phase for E-polarization (top) and H-polarization (bottom).
24
W. Sarakorn / Journal of Applied Geophysics 139 (2017) 16–24
in terms of CPU time are slower than those of the FE method with an unstructured triangular mesh because the processes to evaluate the elemental matrix and coefficient matrix are more complicated. When evaluating the reliability of the codes, it was found that our responses and those calculated by the COMMEMI projects (Zhdanov et al., 1997) are comparable for both polarizations. This shows that our codes are accurate and reliable. Furthermore, our FE codes using an unstructured quadrilateral mesh provide additional power features including local mesh refinement and the capability of handling complex structures, which are useful for a complex 2-D domain. This indicates that our FE codes are flexible and appropriate for use in any real 2-D domain. In the future, these forward modeling codes can be used as a forward calculator to further develop 2-D magnetotelluric inversion for interpreting real magnetotelluric data. Acknowledgments This work was supported by the DPST Research grant [grant number 002/2557].
References Bern, M., 1997. Quadrilateral meshing by circle packing. Int. J. Comput. Geom. Appl. 7–20. Blacker, T.D., Stephenson, M.B., 1991. Paving. A new approach to automated quadrilateral mesh generation. Int. J. Numer. Methods Eng. 32 (4), 811–847.
Franke, A., Börner, R.U., Spitzer, K., 2007. Adaptive unstructured grid finite element simulation of two-dimensional magnetotelluric fields for arbitrary surface and seafloor topography. Geophys. J. Int. 171 (1), 71–86. Key, K., Weiss, C., 2006. Adaptive finite-element modeling using unstructured grids: the 2d magnetotelluric example. Geophysics 71 (6), G291–G299. Lee, S.K., Kim, H., Song, Y., Lee, C.K., 2009. MT2DInvmatlab-a program in MATLAB and FORTRAN for two-dimensional magnetotelluric inversion. Comput. Geosci. 35 (8), 1722–1734. Malleswari, D., Veeraswamy, K., 2014. Numerical simulation of coast effect on magnetotelluric measurements. Acta Geodaet. Geophys. 49 (1), 17–35. Owen, S.J., 1998. A Survey of Unstructured Mesh Generation Technology. International Meshing Roundtable pp. 239–267 Pandey, D., Sinha, M., MacGregor, L., Singh, S., 2008. Ocean coast effect on magnetotelluric data: a case study from Kachchh, India. Mar. Geophys. Res. 29 (3), 185–193. Rao, K.P., Babu, G.A., 2006. EMOD2d–a program in C++ for finite difference modelling of magnetotelluric TM mode responses over 2D earth. Comput. Geosci. 32 (9), 1499–1511. Reddy, I.K., Rankin, D., 1975. Magnetotelluric response of laterally inhomogeneous and anisotropic media. Geophysics 40 (6), 1035–1045. Remacle, J.F., Lambrechts, J., Seny, B., Marchandise, E., Johnen, A., Geuzainet, C., 2012. Blossom-quad: a non-uniform quadrilateral mesh generator using a minimum-cost perfect-matching algorithm. Int. J. Numer. Methods Eng. 89 (9), 1102–1119. Sarrate, J., Huerta, A., 2000. Efficient unstructured quadrilateral mesh generation. Int. J. Numer. Methods Eng. 49 (10), 1327–1350. Singh, U.K., Kant, Y., Singh, R.P., 1995. Effect of coast on magnetotelluric measurements in India. Ann. Geophys. 38 (3-4). Wannamaker, P.E., Stodt, J.A., Rijo, L., 1986. Two-dimensional topographic responses in magnetotellurics modeled using finite elements. Geophysics 51 (11), 2131–2144. Zhdanov, M.S., Varentsov, I.M., Weaver, J.T., Golubev, N.G., Krylov, V.A., 1997. Methods for modelling electromagnetic fields results from COMMEMI - the international project on the comparison of modelling methods for electromagnetic induction. J. Appl. Geophys. 37 (3-4), 133–271.