| Citation: | Yixin Ye, Xiangyun Hu, Dong Xu. A goal-oriented adaptive finite element method for 3D resistivity modeling using dual-error weighting approach. Journal of Earth Science, 2015, 26(6): 821-826. doi: 10.1007/s12583-015-0598-8
|
Finite element (FE) method is a key tool for direct current (DC) resistivity modeling and inversion. A number of FE numerical solutions for 3D DC resistivity problem have been discussed over the past decade (Blome et al., 2009; Huang et al., 2003; Ruan et al., 2001; Sasaki, 1994; Pridmore et al., 1981; Coggon, 1971). Recently, several important uses of FE method such as singularity removal techniques (Li and Spitzer, 2002; Xu et al., 1994), rapid solver (Wu and Wang, 2003; Zhou and Greenhalgh, 2001), anisotropy (Wang et al., 2013; Li and Spitzer, 2005), and surface topography (Penz et al., 2013; Qiang and Luo, 2007; Rücker et al., 2006) have been discussed.
However, for 3D DC resistivity modeling with FE method, regular mesh is commonly employed (Wu and Wang, 2003; Li and Spitzer, 2002) may not accurately accommodate the complex structures such as dipping interfaces and topography. Moreover, model domains are often subdivided according to user's experiences. For simple structures, one can design an optimal mesh by experience, but for complex structures, the design of a reasonable mesh may present difficulty. Due to this, more attention has been paid to the adaptive FE method (Tang et al., 2011; Ren and Tang, 2010). It starts with a coarse mesh and successively refined meshes are generated according to a posteriori error estimator, which assesses the accuracy of FE solution. In the literature, the notion of a posteriori error estimator and adaptivity were first introduced into FE community in 1978 by Babuška and Rheinboldt. Since late 1980's, many methods have been developed including explicit/implicit error estimators, residual type error estimators, recovery type estimators, and error estimators on hierarchic basis. Ovall(2005, 2004) developed a goal-oriented posteriori error estimator based on gradient recovery operator using dual-error weighting approach (DEW). The dual-error weighting approach applies a weighting term to a local error indicator, where the weight is determined by a dual solution of the FE system. It has been applied to the 2D magnetotelluric (MT) problem (Li and Pek, 2008; Key and Weiss, 2006) and the 2.5D controlled-source electromagnetic (CSEM) problem (Li and Key, 2007).
In this paper, we present an adaptive finite element approach for 3D DC resistivity modeling. In the approach, the model domain is subdivided using unstructured tetrahedral mesh which readily allows us to complex structures. Moreover, a goal-oriented posteriori error estimator using dual-error weighting approach is employed to make the mesh more reasonable. After that, we validate the adaptive FE approach on two synthetic resistivity models with analytical solutions and available results from integral equation method. Finally, the applicability of the adaptive FE approach is illustrated on a resistivity model with a topographic ridge.
We assume a 3D conductivity model σ(x, y, z) in a Cartesian system of coordinates (x, y, z) with the origin at the air-earth interface and z positive downwards. The boundary value problem of the 3D DC electrical potential is given by (Xu, 1994)
| \begin{array}{l} \nabla \cdot(\sigma \nabla u)=-I \delta\left(x-x_{s}\right) \delta\left(y-y_{s}\right) \delta\left(z-z_{s}\right) \quad \in \Omega \\ \frac{\partial u}{\partial {\mathit{\boldsymbol{n}}}}=0 \quad \in \Gamma_{s} \\ \frac{\partial u}{\partial {\mathit{\boldsymbol{n}}}}+\frac{\cos (\boldsymbol{r}, \boldsymbol{n})}{|{\mathit{\boldsymbol{r}}}|} u=0 \quad \in \Gamma_{\infty} \end{array} | (1) |
where Ω denotes the model domain, Γs and Γ∞ denote the air-earth interface and the distant boundary respectively, n is the outward normal vector, r is the vector from the source point to the distant boundary (Γ∞) of the model domain Ω, I is the current source located at (xs, ys, zs), and δ is the Dirac delta function.
The use of the weighted residual method (Zienkiewicz and Taylor, 2000) leads to the weak formulation of the boundary value problem
| \begin{array}{l} \int_{\Omega} \sigma \nabla u \cdot \nabla \delta u \mathrm{~d} \Omega+\int_{\Gamma_{\infty}} \sigma \frac{\cos (\boldsymbol{r}, \boldsymbol{n})}{|\boldsymbol{r}|} u \cdot \delta u \mathrm{~d} \Gamma= \\ -\int_{\Omega} I \delta\left(x-x_{s}\right) \delta\left(y-y_{s}\right) \delta\left(z-z_{s}\right) \delta u \mathrm{~d} \Omega \end{array} | (2) |
where δu is an arbitrary variation. For finite element method, the model domain is subdivided into tetrahedral elements. We assume that in each tetrahedral element the electric potential varies linearly and can be approximated by a linear interpolation function. Then the integrals over all the tetrahedron elements can be assembled together to form the global FE system of equations
| {\mathit{\boldsymbol{Ku=P}}} | (3) |
where K is a nd×nd symmetric matrix, nd is the total number of nodes, u is a vector of the unknown potentials at all nodes, P is a nd dimensional vector related to the function Iδ(x-xs)δ(y-ys)δ(z-zs). Solving Eq. (3), the potentials at all nodes are obtained.
The compressed row storage scheme (Barrett et al., 2006) is adopted to deal with the sparse matrix K. For solving the linear system arising from the FE equations, we utilize the matrix-free quasi-minimal residual (QMR) approach (Weiss, 2001), which does not require explicit storage of the coefficient matrix and, hence, requires less memory storage. After the electrical potential u is calculated, apparent resistivity for different configurations can be obtained readily. For the details of tetrahedron element integration, please refer to Rücker et al. (2006).
An adaptive finite element method starts with a coarse mesh, and successively refined meshes are automatically generated according to a posteriori error estimator, which assesses the accuracy of the FE solution. In this paper, we use the dual-error weighting approach (Ovall, 2005, 2004) in which the global influence on the local error is taken into consideration. This method has been recently shown to be effective for the 2D MT (Key and Weiss, 2006) and marine CSEM forward problem (Li and Key, 2007).
First, a coarse mesh is generated by using TetGen (Si, 2003), which tries to force all radius-to-edge ratio below a certain quality constraint. We choose a radius-to-edge of 1.4 for all meshes used for our model tests. After the solution on the coarse mesh is obtained, we use the L2-norm of the difference between a superconvergent gradient recovery operator R=SmQh presented by Bank and Xu (2003) and the piecewise constant gradient to describe a local error indicator for a given element e.
| \eta_{e}=\left\|\left(S^{m} Q_{h}-\mathrm{I}\right) \nabla u_{h}\right\|_{L_{2}(e)} | (4) |
where I is the identification operator, \nabla u_his the FE solution gradient, Qh is the L2 projection operator, S is an appropriate smoothing operator, and m is appropriate number of iterations of the smoothing process, which is chosen to be two for our model tests.
The grid refinement based on the local error indicator (Eq. 4) is suitable for generating a globally accurate solution, hence we will refer to this method as the global refinement. In many cases, however, the solution accuracy is required only at a few selected locations (i.e., the receivers) within the model domain, and several refinement iterations may have little impact on the solution accuracy at these possibly distant locations, as shown in Key and Weiss (2006). In another way, manually increasing the refinement in the proximity of the receivers (i.e., local refinement) would not necessarily work because the solution for the elliptic problem depends on data throughout the entire domain. An efficient solution to this dilemma is to apply a weighting term to the error indicator, where the weight is determined by a dual solution of the FE system, referred to as dual-error weighting approach (Key and Weiss, 2006; Ovall, 2005, 2004).
Consider a functional G that measures the solution error u-uh, where u is the true solution of the partial differential equation and uh is the finite element approximation. Since u is generally unknown, the dual-error weighting approach uses a dual problem to approximate G. The dual problem can be expressed as
| B^{*}(w, v)=G(v) | (5) |
where the operator B denotes the operations on the left side of Eq. (2), B* is a dual operator and is defined as B*(w, v)=B(w, v), and w=u, v=δu.
To compute the FE solution wh of the dual problem, we need to solve the dual problem with a sensible choice for the error functional G. We use the seminorm of uh to approximate G (Ovall, 2004). The dual problem is then written as
| B^{*}(w, v)=G(v) \equiv \int_{\Omega_{s}}(R-I) \nabla u_{h} \cdot \nabla v \mathrm{~d} \Omega_{s} | (6) |
After the FE solution wh is computed, the dual-error weighting indicator for each element is defined as
| \hat{\eta}_{e}=\eta_{e} \bar{\eta}_{e} | (7) |
with \bar{\eta}_{e}=\sigma\left\|(R-I) \nabla w_{h}\right\|_{L_{2}(e)}.
Finally, the mesh is refined adaptively according to the \hat η_e value. In our implementation, we found that refining 10% of elements with the largest \hat η_e values gives a good rate of convergence without over-refining the grid or greatly increasing the computational time per iteration. After an improved mesh is generated, the problem is solved on the new mesh and the error indicator \hat η_e is updated.
In the modeling, the mesh refining process repeats iteratively until the maximum relative difference between the current and previous grid's FE solution at all the receiver locations is less than 1%. All the model computations are performed on the personal computer with double-cored 2.9 GHz CPU and 8 GRAM.
First, we consider the case of two vertical layered model (Fig. 1a), with dimensions of 2 000×2 000×1 000 m. The left vertical layer has a resistivity of 10 Ω·m, the right vertical layer has a resistivity of 1 000 Ω·m. To discuss the accuracy of numerical solutions, the pole-pole sounding along the y-axis is carried out with the single source electrode A located at -10 m on the surface. Measurements are taken by 28 potential electrodes ranging from -9 to 18 m along the y-axis with 1 m electrode spacing. The initial mesh generated by the TetGen consists of 679 nodes and 3 024 elements (Fig. 1b).
As the algorithm is performed, the initial mesh refined to an intermediate mesh which consists of 3 004 nodes and 15 971 elements after four iterations (Fig. 1c), it then refined to a final mesh which consists of 9 102 nodes and 51 369 elements after seven iterations (Fig. 1d). From Fig. 1, we can see the mesh is automatically refined in the area of potential electrodes especially in the vicinity of the source electrode, because the element error percentage indicates their singularity, and the adaptive strategy refines the corresponding regions. The refinement iterations took about 0.13 min, while the total computation took about 0.37 min. We compare the numerical solutions for three meshes with the analytical solutions (Figs. 2 and 3). Figure 2 shows that the numerical solutions can asymptotically converge to analytical solutions with mesh refinements. Figure 3 shows that the accuracy of numerical solutions is dramatically improved with mesh automatically refined especially in the vicinity of the source electrode. On the initial mesh, in the vicinity of the source electrode we observe quite large error 16.5%, whereas with increasing distance the errors stabilize at around 4%. On the final mesh, an average relative deviation of 0.5% and 1.8% near the source electrode is obtained.
The second model is a cube buried in a two vertical layered earth (Fig. 4a). The left vertical layer has a resistivity of 10 Ω·m, the right vertical layer has a resistivity of 100 Ω·m. A conductive cube with a resistivity of 10 Ω·m and a side length of 2 m is buried near the contact interface. A Schlumberger sounding is carried out along the y-axis and the two current sources are located at (0, -3.4 m, 0) and (0, -2.6 m, 0). The apparent resistivity is calculated with half length of potential electrodes (AB/2) varying from 1 to 100 m.
The initial mesh generated by the TetGen consists of 775 nodes and 3 525 elements (Fig. 4b). As the algorithm is performed, the mesh is automatically refined, which results in a intermediate mesh of 4 751 nodes and 26 275 elements after six iterations (Fig. 4c) and a final mesh of 54 666 nodes and 330 430 elements after twelve iterations (Fig. 4d). From Fig. 4, we can see that the areas near each receiver location especially in the vicinity of the source electrodes have been automatically refined. The refinement iterations took about 2.2 min, while the total computation took about 6.5 min. The numerical apparent resistivity and relative deviation curves in comparison with the integral equation results from a boundary integral method (Hvoždara and Kaikkonen, 1994) are shown in Figs. 5 and 6. Figure 5 demonstrates that the numerical apparent resistivity values in the adaptive process can asymptotically converge to the integral equation results, and Fig. 6 shows the relative deviations decrease largely in the adaptive process. On the initial mesh, the maximal relative deviation is about 10% and the average relative deviation is 4%, while on the final mesh, we get rather accurate results with maximal relative deviation of 1.2% and an average relative deviation of 0.4%.
The last model is a 3D topographic ridge with a resistivity block embedded, with its surface view and cross section shown in Fig. 7. The strike of the ridge is 20 m along the x-axis. The half space has a resistivity of 20 Ω·m, with dimensions of 2 000×2 000×1 000 m. A resistive block of 500 Ω·m with dimensions of 20, 4 and 30 m along x-, y- and z-axes, respectively. A mid-gradient configuration perpendicular to the strike is carried out with double source electrodes A and B located at -80 and 80 m respectively. Measurements are taken by 30 electrodes located from -29 to 29 m with 2 m spacing of each electrodes.
Figure 8 shows the cross section and the surface view of final refined mesh after twelve iterations, which contains 50 396 nodes and 302 687 elements. It shows that the source electrodes area and the topographical regions are effectively refined. The refinement iterations took about 2 min, while the total computation took about 6 min. On the final refined mesh, the total apparent resistivity response (ρs) and the pure topographic apparent resistivity response (ρsD) are calculated sepa-rately, then the apparent resistivity response with topographic correction (ρsD) is estimated using the ratio method
| \rho_{s}^{G}=\rho_{s} / \frac{\rho_{s}^{D}}{\rho_{1}} | (8) |
where ρ1 is the resistivity of the half space. It is then compared with the apparent resistivity with flat surface.
Figure 9 displays the different apparent resistivity curves of the 3D topographic ridge using the mid-gradient configuration. It shows that the total apparent resistivity curve is similar to the pure topographic apparent resistivity curve, this means that the total apparent resistivity curve is strongly influenced by the topographic ridge. And the apparent resistivity curve with topographic correction is similar to the apparent resistivity curve with flat surface, which shows the anomaly of high resistivity block correctly.
This paper presents a goal-oriented adaptive FE method to the 3D DC modeling problem. Using the goal-oriented concept, the accuracy of potentials at electrodes can be significantly improved. The validations show that the adaptive FE algorithm can offer a reliable numerical solution, which can asymptotically converge to accurate solutions with grid refinements. The unstructured tetrahedral mesh can easily accommodate complex structures, such as dipping interfaces and rough topography. It is seen from the visualized meshes, the areas near each receiver location especially in the vicinity of the source electrodes have been automatically refined. The modeling results of models 1 and 2 illustrate the kinds of structures that can be accurately modeled with this code. Only about 30% of the total computational time was spent on the refinement iterations. Model 3 study shows a clearly non-negligible effect of the topography in apparent resistivity data. More powerful properties of the adaptive refinement methodology could be expected in geometrically complex cases.
ACKNOWLEDGMENTS: This study was supported by the National Natural Science Foundation of China (No. 41204055), the National Basic Research Program of China (No. 2013CB733203), and the Opening Project (No. SMIL-2014-06) of Hubei Subsurface Multi-Scale Imaging Lab (SMIL), China University of Geosciences, Wuhan, China. The authors give special thanks to Prof. Yuguo Li for help in preparing materials and critical suggestions on the manuscript. We thank Hang Si for the grid generation and refinement code TetGen, Mi Hvozdara for integral equation code, and Chester Weiss for the matrix-free QMR code utilized in this paper. We would also thank Prof. Changchun Yin and two anonymous reviewers for their helpful comments and constructive suggestions.| Babuška, I., Rheinboldt, W. C., 1978. A-Posteriori Error Estimates for the Finite Element Method. International Journal for Numerical Methods in Engineering, 12(10): 1597–1615. doi: 10.1002/nme.1620121010 |
| Bank, R. E., Xu, J. C., 2003. Asymptotically Exact a Posteriori Error Estimators, Part Ⅱ: General Unstructured Grids. SIAM Journal on Numerical Analysis, 41: 2313–2332 doi: 10.1137/S0036142901398751 |
| Barrett, R., Berry, M., Chan, T. F., et al., 2006. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Society for Industrial and Applied Mathematics, Philadelphia |
| Blome, M., Maurer, H. R., Schmidt, K., 2009. Advances in Three-Dimensional Geoelectric Forward Solver Techniques. Geophysical Journal International, 176(3): 740–752. doi: 10.1111/j.1365-246x.2008.04006.x |
| Coggon, J. H., 1971. Electromagnetic and Electrical Modeling by the Finite Element Method. Geophysics, 36(1): 132–155 doi: 10.1190/1.1440151 |
| Huang, J. G., Ruan, B. Y., Bao, G. S., 2003. Finite Element Method for IP Modeling on 3-D Geoelectric Section. Earth Science—Journal of China University of Geosciences, 28(3): 323–326 (in Chinese with English Abstract) http://en.cnki.com.cn/Article_en/CJFDTOTAL-DQKX200303013.htm |
| Hvoždara, M., Kaikkonen, P., 1994. The Boundary Integral Calculations of the Forward Problem for D.C. Sounding and MMR Methods for a 3-D Body near a Vertical Contact. Studia Geophysica et Geodætica, 38(4): 375–398. doi: 10.1007/bf02296169 |
| Key, K., Weiss, C., 2006. Adaptive Finite-Element Modeling Using Unstructured Grids: The 2D Magnetotelluric Example. Geophysics, 71(6): G291–G299. doi: 10.1190/1.2348091 |
| Li, Y. G., Key, K., 2007. 2D Marine Controlled-Source Electromagnetic Modeling: Part 1—An Adaptive Finite-Element Algorithm. Geophysics, 72(2): WA51–WA62. doi: 10.1190/1.2432262 |
| Li, Y. G., Pek, J., 2008. Adaptive Finite Element Modelling of Two-Dimensional Magnetotelluric Fields in General Anisotropic Media. Geophysical Journal International, 175(3): 942–954. doi: 10.1111/j.1365-246x.2008.03955.x |
| Li, Y. G., Spitzer, K., 2005. Finite Element Resistivity Modelling for Three-Dimensional Structures with Arbitrary Anisotropy. Physics of the Earth and Planetary Interiors, 150(1–3): 15–27. doi: 10.1016/j.pepi.2004.08.014 |
| Li, Y. G., Spitzer, K., 2002. Three-Dimensional DC Resistivity Forward Modelling Using Finite Elements in Comparison with Finite-Difference Solutions. Geophysical Journal International, 151(3): 924–934. doi: 10.1046/j.1365-246x.2002.01819.x |
| Ovall, J. S., 2004. Duality-Based Adaptive Refinement for Elliptic: [Dissertation]. University of California, San Diego |
| Ovall, J. S., 2005. Asymptotically Exact Functional Error Estimators Based on Superconvergent Gradient Recovery. Numerical Mathematics, 102: 543–558 |
| Penz, S., Chauris, H., Donno, D., et al., 2013. Resistivity Modelling with Topography. Geophysical Journal International, 194(3): 1486–1497 doi: 10.1093/gji/ggt169 |
| Pridmore, D., Hohmann, G. W., Ward, S. H., et al., 1981. An Investigation of Finite-Element Modeling for Electrical and Electromagnetic Modelling Data in Three Dimensions. Geophysics, 46: 1009–1024 doi: 10.1190/1.1441239 |
| Qiang, J. K., Luo, Y. Z., 2007. The Resistivity FEM Numerical Modeling on 3D Undulating Topography. Chinese J. Geophys. , 50(5): 1606–1613 (in Chinese with English Abstract) |
| Ren, Z. Y., Tang, J. T., 2010. 3D Direct Current Resistivity Modeling with Unstructured Mesh by Adaptive Finite-Element Method. Geophysics, 75(1): H7–H17 doi: 10.1190/1.3298690 |
| Ruan, B. Y., Xiong, B., Xu, S. Z., 2001. Finite Element Method for Modeling Resistivity Sounding on 3D Geoelectric Section. Earth Science—Journal of China University of Geosciences, 26(1): 73–77 (in Chinese with English Abstract) |
| Rücker, C., Günther, T., Spitzer, K., 2006. Three-Dimensional Modelling and Inversion of DCResistivity Data Incorporating Topography-I. Modelling. Geophysical Journal International, 166(2): 495–505. doi: 10.1111/j.1365-246x.2006.03010.x |
| Sasaki, Y., 1994. 3-D Resistivity Inversion Using the Finite-lement Method. Geophysics, 59(12): 1839–1848. doi: 10.1190/1.1443571 |
| Si, H., 2003. TETGEN: A 3D Delaunay Tetrahedral Mesh Generator. http://tetgen.berlios.de. [2015.11.10] |
| Tang, J. T., Wang, F. Y., Xiao, X., et al., 2011. 2.5-D DCResistivity Modeling Considering Flexibility and Accuracy. Journal of Earth Science, 22(1): 124–130 doi: 10.1007/s12583-011-0163-z |
| Wang, W., Wu, X. P., Spitzer, K., 2013. Three-Dimensional DC Anisotropic Resistivity Modelling Using Finite-Elements on Unstructured Grids. Geophysical Journal International, 193(2): 734–746. doi: 10.1093/gji/ggs124 |
| Weiss, C. J., 2001. A Matrix-Free Approach to Solving the Fully 3D Electromagnetic Induction Problem. 71st Annual International Meeting, SEG, Expanded Abstracts |
| Wu, X. P., Wang, T. T., 2003. A 3-D Finite Element Resistivity Forward Modeling Using Conjugate Gradient Algorithm. Chinese J. Geophys. , 46(3): 428–432 (in Chinese with English Abstract) |
| Xu, S. Z., 1994. The Finite Element Method in Geophysics. Science Press, Beijing (in Chinese) |
| Xu, S. Z., Liu, B., Ruan, B. Y., 1994. The Finite Element Method for Solving Anomalous Potential for Resistivity Surveys. Chinese J. Geophys. , 37(S2): 511–515 (in Chinese with English Abstract) |
| Zhou, B., Greenhalgh, S. A., 2001. Finite Element Three-Dimensional Direct Current Resistivity Modelling: Accuracy and Efficiency Considerations. Geophysical Journal International, 145: 679–688 doi: 10.1046/j.0956-540x.2001.01412.x |
| Zienkiewicz, O. C., Taylor, R. L., 2000. The Finite Element Method (5th Ed. ), Basic Foundation. Butterworth-Heinemann, Elsevier |
| 1. | Lichao Nie, Zhaoyang Deng, Zhi-Qiang Li, et al. Detection method for water-bearing structure in tunnel by DC resistivity based on directional indicating variable mesh. Tunnelling and Underground Space Technology, 2024, 149: 105820. doi:10.1016/j.tust.2024.105820 | |
| 2. | Supti Sadhukhan. 3D Pore Geometry and Electrical Conductivity. Geotechnical and Geological Engineering, 2024. doi:10.1007/s10706-024-02803-7 | |
| 3. | Yi’an Cui, Zhixue Chen, Xiaoxiong Zhu, et al. Sequential and simultaneous joint inversion of resistivity and IP sounding data using particle swarm optimization. Journal of Earth Science, 2017, 28(4): 709. doi:10.1007/s12583-017-0749-1 | |
| 4. | Jianchao Cai, Wei Wei, Xiangyun Hu, et al. Electrical conductivity models in saturated porous media: A review. Earth-Science Reviews, 2017, 171: 419. doi:10.1016/j.earscirev.2017.06.013 |
Figures(9)
Copyright © 2013-2020 Journal of Earth Science 鄂ICP备15021562号-2
Tel: +86-27-67885075 Fax: +86-27-67885075 E-mail: xbb@cug.edu.cn
Address: Editorial Office of Journal, China University of Geosciences, Yujiashan, Wuhan, Hubei 430074, P. R. China
Supported by:
Beijing Renhe Information Technology Co. Ltd
E-mail:
info@rhhz.net
Yixin Ye, Xiangyun Hu, Dong Xu. A goal-oriented adaptive finite element method for 3D resistivity modeling using dual-error weighting approach. Journal of Earth Science, 2015, 26(6): 821-826. doi: 10.1007/s12583-015-0598-8
DownLoad:
DownLoad:
