TY - JOUR
AU - Zhi,, Qingquan
AB - Abstract Traditional nuclear magnetic resonance (NMR) is mainly applied to horizontal ground and one-dimensional electrical structures. However, calculations of the excitation fields rarely consider the three-dimensional electrical changes of the subsurface medium and undulating terrain, or the deformations caused by emission sources. Therefore, to analyse the influences of terrain fluctuations, emission source deformations and three-dimensional electrical changes on NMR, three-dimensional finite element forward modelling of undulating terrain NMR was conducted in this study. First, based on a scalar finite element method, the direct calculations of the excited magnetic fields of a three-dimensional electrical medium were realised, which improved calculation accuracy by avoiding the finite element calculations of magnetic vector potential and vector on the magnetic field. During the source loading process, the equivalent thin wire source of the pseudo δ function was used to load the source function directly into the equation for the purpose of achieving total field calculations. This was completed to enable the calculations to be applied to any shape of the transmitting loop and undulating terrain. Then, the components of the excitation magnetic fields perpendicular to the geomagnetic fields were calculated using the rotation matrix. Finally, the NMR sensitivity function and 3D responses were calculated. The calculations of the excitation magnetic fields were verified using a uniform half-space model. The overall algorithm was tested by the nuclear magnetic responses of the layered medium. Also, a typical undulating terrain model was adopted and the complex excitation source NMR was simulated using the algorithm proposed in this study. The algorithm provided a three-dimensional forward basis for the NMR inversion in the cases of determining the electrical medium for the subsequent undulating terrain. nuclear magnetic resonance, three-dimensional finite element, underground electromagnetic field, complex field source morphology, inhomogeneous medium 1. Introduction The ground nuclear magnetic resonance (NMR) method is a geophysical method that involves laying ground coils for the excitation of low-frequency electromagnetic fields (excitation frequency about ∼2–3 kHz) to measure the spin echo responses of underground hydrogen atoms and search directly for water sources. This method is capable of determining the depth, thickness and content of the aquifer during an underground water source search process, and has been proven to be effective. In the traditional NMR forward simulation, the calculation of the excitation fields usually adopts the analytical solutions of the horizontal ground uniform electric medium electromagnetic fields (Lin 2010), or the integral filtering solutions of horizontal layered medium (Wong et al. 2007). These methods are generally only applicable to horizontal terrain, as well as regular transmitting coils and one-dimensional resistivity structures. Therefore, to avoid these limitations and improve the applicability of the forward algorithms to arbitrary topography, irregular transmitting coils and three-dimensional electrical media, a three-dimensional finite element forward algorithm for NMR was proposed in the current study. This algorithm was based on the direct calculations of the excitation magnetic fields. It has been found that the three-dimensional simulations of active Maxwell's Equations are difficult to perform. Previously, Ryu et al. (1970) provided the integral expressions of electromagnetic field components generated by loop sources in a horizontally layered medium. Since this time, active electromagnetic field calculations have consistently been topics of interest in the field of geophysics. In particular, controllable-source electromagnetic detections, such as control-source magnetotelluric (CSMT) and transient electromagnetics (TEM) have been thoroughly examined. During the calculations of the influences of emission loop shapes on electromagnetic fields, Raiche (1987) had presented the results of the electromagnetic responses of a polygonal transient electromagnetic layered medium. Also, Cui & Chew (2000) proposed a forward modelling method of liner antennae with arbitrary shapes and directions on layered ground. Zhou et al. (2011) used a dipole superposition principle to realise the calculations of the electromagnetic responses of arbitrary polygon loop sources. Li & Han (2017) discussed the influences of complex field sources on the electromagnetic responses of ocean controllable sources based on a staggered-grid finite volume method. Li et al. (2018) realised the transient electromagnetic simulations of complex morphological loops using a time domain finite element method. However, due to the changes in the shapes of the emission sources, the distributions of underground electromagnetic fields will inevitably be affected. Subsequently, the distributions of the magnetic resonance kernel functions will also be affected, and the magnetic resonance responses will ultimately become impacted. Therefore, there is a need to further investigate the three-dimensional finite element forward modelling of complex field source shapes and undulating surface NMR. The three-dimensional forward numerical simulation of electromagnetic fields includes an integral equation method, finite difference method, finite volume method and finite element method. Previously, Wannamaker et al. (1984) realised the electromagnetic response calculations of three-dimensional anomalous bodies in the layered medium backgrounds based on an integral equation. Xiong (1992) achieved the iterative integral forward modelling of electromagnetic fields with three-dimensional geo-electric structures, which had effectively reduced the memory requirements during the forward process. Zhdanov et al. (2007) presented an integral equation method for the 3D electromagnetic modelling of complex structures with non-uniform background conductivity. In another related study. Avdeev & Knizhnik (2009) improved the algorithm based on an integral equation, and then verified the algorithm using an induction logging model. It has been found that for NMR forward modelling it is necessary to accurately calculate the magnetic fields of each grid node in an underground site. This method differs from the traditional controllable-source electromagnetic method and the transient electromagnetic method, since those methods only require the calculations of the electromagnetic fields at the measuring point positions. Therefore, it has been found that integral equation methods are unable to properly apply the requirements of surface NMR. In comparison, such global singulation methods as finite difference, finite volume and finite element methods have been found to be more suitable for the requirements of NMR forward modelling. For example, Newman & Alumbaugh (1995) adopted a finite difference method for the forward modelling and sensitivity matrix calculations during the realisation processes of controllable-source three-dimensional inversions. Mogi (1996) realised the three-dimensional forward modelling of magnetotelluric electromagnetics based on a staggered grid finite difference method. In the studies conducted by Jaysaval et al. (2014), a fast finite difference forward modelling method of multi-model 3D electromagnetic field was successfully achieved. Haber & Ascher (2001) and Haber et al. (2007) proposed an effective algorithm for realizing the three-dimensional forward modelling of electromagnetic fields using a finite volume method. Also, Jahandari & Farquharson (2014) realized the solution of the electromagnetic forwarding problem by using an unstructured grid finite volume method. In the previously mentioned studies, it had been found that as long as the local variations in conductivity and surface topography were not very severe, the finite difference and finite volume methods had worked well. However, if the surface undulations and electrical structure changes were irregular, it was observed that the finite difference method and the finite volume method required additional complex processing techniques or increases in meshing fineness, which had resulted in the algorithm implementations becoming relatively complex. It has been found that finite element methods are more suitable for the calculations of electromagnetic fields in strong non-uniform medium structures and undulating surfaces. In the studies conducted by Haber (1999) and Haber et al. (2007) the three-dimensional forward modelling of frequency domain electromagnetic fields was successfully realized based on a finite element method. Also, Li et al. (2011, 2013) calculated the anomaly fields of frequency-domain electromagnetic fields based on a vector finite element method, and Cai et al. (2017a, b) calculated the total field of controllable source magnetotellurics based on a finite element method. A time domain vector finite element method was adopted by Li et al. (2018) for the purpose of obtaining the three-dimensional forward modelling algorithms of the transient electromagnetic fields of complex shaped emission wireframes. In previous studies, to avoid pseudo-solution problems, the three-dimensional finite element forward modelling of active electromagnetic fields have mainly adopted vector finite element methods. Also, for the purpose of avoiding the discontinuity influences of the resistivity, the electric field components are often calculated first, and then the magnetic fields are calculated using Maxwell's Equations. During the calculations of magnetic fields, the introduction of differentials often results in low accuracy in the magnetic field calculations. Lehmann-Horn et al. (2011) presented a hybrid finite element algorithm to realise the three-dimensional forward modelling of NMR. First, an anomalous field algorithm was adopted to realise the calculations of the electromagnetic fields of arbitrary emission shape loop sources using line integration, and the magnetic vector potential was calculated using scalar finite element methods. The total magnetic fields of subsurface grid nodes can be calculated from the magnetic vector potential. However, the calculations of magnetic fields using the magnetic vector potential also requires the solution of the magnetic fields, which may potentially lead to lower accuracy rates of the magnetic field calculations. In addition, the calculations of the primary fields using the line integrals are often limited by the terrain, and it may not be possible to calculate the analytical solutions of arbitrary undulating terrain. Although previous studies have listed cases involving undulating terrain, the undulating terrain in those cases were equivalent to the overlaying of two horizontal terrains. Therefore, to avoid the derivation of the electrical fields and magnetic vector potential in the calculations of magnetic fields for the solutions of vector finite elements, as well as the complex limitations of primary field calculations in abnormal field algorithms, a three-dimensional finite element forward calculation of direct magnetic fields based on the total field was proposed in this study. The excitation sources were directly added into the equation, and the scalar finite elements were used to solve the variational problems with divergence constraint. The total magnetic field was directly calculated, which avoided the error introduced by the derivation and successfully avoided the limitations of calculating the primary fields. It has been determined that NMR signals are weak and susceptible to interference. Shushakov (1996) pointed out that in the cases of large background noise, the signal-to-noise ratios can be improved to some extent by using ‘8’-shaped coils. To reduce the interference and improve the signal-to-noise ratio, some previous studies have implemented separate coil technology. However, there have been few analyses completed regarding the effects of different emission device parameters on NMR. In addition, the influences of NMR on the distributions of electrical media have recently attracted the attention of many researchers (Legchenko et al.2008). For typical problems related to electrical influence, Lin (2010) pointed out that average resistivity methods can be used for simplification purposes when the resistivity of the cover layer changes. In the studies conducted by Wong et al. (2007), the laws of NMR responses under layered medium conditions were analysed and it was pointed out that conductivity was a prerequisite for generating ground NMR phase signals. Braun et al. (2005) and Braun & Yaramanci (2008) attempted to obtain electrical distributions by using the phase information of NMR signals and had achieved good results during theoretical model testing. Hertrich & Yaramanci (2002) and Hertrich et al. (2005) used an electrical sounding method to obtain the electrical model of subsurface medium. Then, based on a simulated annealing method, the joint inversion of a NMR method and an electrical sounding method was realised. Legchenko et al. (2009) and Wan et al. (2013) adopted genetic algorithm for transient electromagnetic and NMR joint inversion. However, the above-mentioned forward and inversion interpretations were all based on one-dimensional layered models as the research targets, without consideration given to such factors as the three-dimensional variations of the dielectric properties and the fluctuations of the terrain. In this study, a three-dimensional finite element forward algorithm was proposed for NMR, which had the ability to simulate the electromagnetic fields of undulating surfaces, complex excitation field sources and three-dimensional geoelectric structures. This problem was represented by a magnetic field. This study's magnetic field equation was established by Maxwell's Equation, and the divergence conditions were imposed by penalty terms to avoid pseudo-solutions. The current-carrying wires were added to the equation using the pseudo-δ source to realise the finite element calculations of the electromagnetic fields that were suitable for arbitrary emission shapes and undulating terrain. The magnetic field distributions were determined to calculate the sensitivity kernel functions of the water-body distributions, thereby solving the complete surface MRI forward modelling problems. A layered model was used to verify the correctness of the algorithm, and the NMR response characteristics of the different emission wireframe shapes were calculated. The results of the new algorithm, which had been applied to a case of three-dimensional mountain karst flushing synthesis, are presented in this to highlight the effects of the surface morphology, conductivity structures and excitation coil shapes on the NMR signals. 2. Finite element forward principle of the NMR excitation fields Water searches are conducted using NMR excited by feeding a sinusoidal current of Larmor frequency into ground-emitting loop-transmitting coils. After the external magnetic fields are removed, a free decay signal equal to the Larmor frequency is released. Then, the initial amplitude of the free decay signal will be dependent on the following equation: $$\begin{equation}{E_0}({\rm{q}}) = \frac{{{\omega _0}{M_0}}}{I}\int{{{B_ \bot }}}\sin \left( {\gamma \frac{{{B_ \bot }}}{I}q} \right)n(r)dv.\end{equation}$$ (1) In the equation, E0 is the initial amplitude of free decay signal; ω0 represents the local Larmor frequency; M0 represents the proton magnetic moment; I is the emission current; γ denotes the proton magnetic rotation ratio; n(r) indicates the water content at the point r; q is the emission pulse moment and B⊥ represents the component of the applied magnetic field perpendicular to the geomagnetic field. The initial amplitude of the free decay signal is related to the moisture content, nuclear magnetic rotation ratio and magnitude of the applied magnetic field along the geomagnetic field component. Therefore, under the premise of determining the geoelectric parameters, to obtain the initial amplitude the calculation of the excitation field B⊥ is required. 2.1. Control equations and boundary conditions By assuming that the time harmonic factor as e−iωt, the active harmonic electromagnetic field can be obtained as follows: $$\begin{equation}\left. {\begin{array}{@{}*{1}{l}@{}} {\nabla \times {{\bf E}} = i\omega \mu {{\bf H}}}\\ {\nabla \times {{\bf H}} = (\sigma - i {\omega} \varepsilon ){\bf {E}} + {{\bf {J}}_s}}\end{array}} \right\}.\end{equation}$$ (2) In the equation, ω represents the circular frequency of the excitation loop source; i is the imaginary unit; t denotes the time; E is the electric field vector; H indicates the magnetic field vector and Js is the current density of the applied current source. For the second side of the above-mentioned formula, the rotation is obtained and simplified and can be written as the following control equation: $$\begin{equation}\nabla \times \nabla \times {\bf {H}} - {k^2}{\bf {H}} = \nabla \times {\bf {J}}.\end{equation}$$ (3) Since the ground nuclear magnetic excitation frequency belongs to the low frequency field, the current can be ignored. Therefore, k2 = iωμσ. When solving equation (3), a node finite element method is often used to calculate the vector potential function. In other cases, the vector finite element is applied to first calculate the electrical field to avoid the calculation of the current density rotation or the discontinuity of the normal electrical field. Then, the magnetic field value is obtained by derivation. However, the magnetic field value calculated by the derivation will have a lower accuracy rate and may even be lower than the accuracy directly solved using a finite element method. Therefore, a finite element method should be used to directly calculate the magnetic field value. In the current study, for simplicity of calculation, the effects of the magnetic field on the outer boundary of infinity were ignored. Then, in accordance with the variation principle, the functional that satisfied equation (3) was obtained as follows: $$\begin{equation}F({\bf {H}}) = \frac{1}{2} \iiint \nolimits_V {\left[ {(\nabla \times {\bf {H}}) \cdot (\nabla \times {\bf {H}}) - {k^2}{\bf {H}} \cdot {\bf {H}}} \right]dV} - \iiint \nolimits_V {{\bf {H}} \cdot (\nabla \times {\bf {J}})dV}.\end{equation}$$ (4) 2.2. Pseudo-solution problems The first-order variation of the functional (4) was taken, which indicated that the test field H should be second-order differentiable. However, in the actual numerical simulation of finite element methods, the discretization requires only that the interpolation function or the expansion function is continuous and does not require its derivative. Therefore, the solution obtained will be a ‘weak solution’ that does not strictly satisfy the governing differential equation. These ‘weak solutions’ are sometimes wrong and do not satisfy the divergence conditions|$\nabla \cdot {\bf {H}} = 0$|⁠. In such cases, divergence conditions are imposed by adding penalty terms, and the functional will then be as follows: $$\begin{equation}F({\bf {H}}) = \frac{1}{2}\iiint\nolimits_V {\left[ {(\nabla \times {\bf {H}}) \cdot (\nabla \times {\bf {H}}) - {k^2}{\bf {H}} \cdot {\bf {H}}} \right]dV} - \iiint\nolimits_V {{\bf {H}} \cdot rot({\bf {J}})dV} + \frac{1}{2}\iiint\nolimits_V {{{(\nabla \cdot {\bf {H}})}^2}dV}.\end{equation}$$ (5) 2.3. Source loading NMR is powered by thin wires with negligible volume, which tends to lead to difficulties in loading sources. In theoretical derivation, for a field source with negligible volume, it is often processed using a function δ. However, in numerical calculations, the field values near field sources are singular. In addition, the calculation of the δ source rotation has been found difficult to achieve. Therefore, traditional methods are used to separately solve the background fields and secondary fields. However, the disadvantages of these methods are that it is often difficult to determine suitable background structures when the surfaces are undulating or the geological bodies are complex. Another simple and effective solution is to use an approximation function to equalize the roles of the field sources. The sources will be distributed within a certain range for the purpose of avoiding singularity at the source points. As a result, the background fields and anomalous fields will not be distinguishable to directly solve the total solution field values. In this study, the pseudo-δ function method was adopted and the function expression is shown in equation (6). In the equations, the parameter τ indicates a parameter that controls the width and amplitude of the source distribution; x is the coordinate of a spatial dimension and x0 represents the location of the source (for example, the coordinate value of a singular point). The effects of thin wire sources can be equivalent to using a pseudo-δ function. For example, the current intensities at the space I(x0, y0, z0) will be the same as the responses of the wire sources along the direction y. Also, the current density distribution will be Jy = Iδs(x − x0)(z − z0). Therefore, by using a pseudo-δ function to describe a source item will be equivalent to a certain degree of broadening of the wire. The effects of this broadening can be ignored at the locations far from the field sources. Therefore, the calculation of the current density in equation (5) can be converted into a derivative operation on the pseudo-δ function as follows: $$\begin{equation}{\delta _s}(x - {x_0}) = \frac{1}{{2\tau }}\left\{ {\begin{array}{@{}*{2}{l}@{}} 0&\quad\quad{(x - {x_0}) \le - 2\tau ,}\\ {[(x - {x_0} + 2\tau )/\tau ]/2}&\quad\quad{ - 2\tau < (x - {x_0}) \le - \tau ,}\\ { - {{[(x - {x_0} + 2\tau )/\tau ]}^2} + 2(x - {x_0} + 2\tau )/\tau - 1}&\quad\quad{ - \tau < (x - {x_0}) \le \tau ,}\\ {[(x - {x_0} + 2\tau )/\tau ]/2 - 4(x - {x_0} + 2\tau )/\tau + 8}&\quad\quad{\tau < (x - {x_0}) \le 2\tau ,}\\ 0&\quad\quad{2\tau < (x - {x_0}).} \end{array}}\right.\end{equation}$$ (6) 2.4. Meshing and unit analysis results In finite element calculations, the calculation areas can be flexibly divided according to the electrical distribution characteristics of the study area using any suitable shape of grid. Then, to adapt to the terrain fluctuations and the three-dimensional distributions of the underground electrical medium, the entire area should be divided by tetrahedral units. Under the premise of ensuring the accuracy of the results, the number of calculation grids was reduced as much as possible. Then, non-uniform meshing can be applied. As shown in figure 1, the loop was placed in the center of the undulating surface. Since it was required to obtain a fine magnetic field value, a finer meshing was adopted. In the boundary regions far from the target, the mesh can be seen to become rapidly sparse in the direction of the outer boundaries to reduce the order of the matrix and minimise the effects of the boundaries. Figure 1. Open in new tabDownload slide Schematic diagram of non-uniform meshing. The meshing schematic diagram shows that the grid has a good fit to the terrain and transmitter loop, and enlarges at the boundary to reduce the amount of calculation. Figure 1. Open in new tabDownload slide Schematic diagram of non-uniform meshing. The meshing schematic diagram shows that the grid has a good fit to the terrain and transmitter loop, and enlarges at the boundary to reduce the amount of calculation. In the current study, the solution area was divided for the purpose of discretising equation (4). Then, by assuming that the internal physical properties of the meshing unit were uniform, the three components of the magnetic field on the element node i were set as Hxi, Hyi, Hzi. The stagnation point for solving the variational problem (5) for all of the elements after meshing was equivalent to solving the following equations: $$\begin{equation}\left\{ {\begin{array}{@{}*{1}{l}@{}} {\frac{{\partial F\left( {{H_x}} \right)}}{{\partial {H_{xi}}}} = \sum\limits_{e = 1}^n {\frac{{\partial {F_e}}}{{\partial H_{xi}^e}}} = 0,}\\ {\frac{{\partial F\left( {{H_y}} \right)}}{{\partial {H_{yi}}}} = \sum\limits_{e = 1}^n {\frac{{\partial {F_e}}}{{\partial H_{yi}^e}}} = 0,}\\ {\frac{{\partial F\left( {{H_z}} \right)}}{{\partial {H_{zi}}}} = \sum\limits_{e = 1}^n {\frac{{\partial {F_e}}}{{\partial H_{zi}^e}}} = 0,} \end{array}} \right.\end{equation}$$ (7) where n is the number of meshing units and e is the unit number, let|$R = \nabla \times {\bf {J}}$|⁠. Then, the linear equations in each unit will be as follows: $$\begin{equation}\left\{ {\begin{array}{@{}*{1}{c}@{}} {\frac{{\partial {F_e}}}{{\partial H_{xi}^e}} = \sum\limits_{j = 1}^{ne} {\int\limits_{e}{{\left[ {\left( {\frac{{\partial N_i^e}}{\partial {y}}\frac{{\partial N_j^e}}{\partial {y}} + \frac{{\partial N_i^e}}{{\partial z}}\frac{{\partial N_j^e}}{{\partial z}} - {k^2}N_i^eN_j^e} \right)H_{xj}^e - \frac{{\partial N_i^e}}{\partial {y}}\frac{{\partial N_j^e}}{{\partial x}}H_{yj}^e - \frac{{\partial N_i^e}}{{\partial z}}\frac{{\partial N_j^e}}{{\partial x}}H_{zj}^e} \right]dxdydz}} = \sum\limits_{j = 1}^{ne} {\int\limits_{e}{{{N_i}{N_j}dxdydz{R_{xj}}}}} ,} }\\ {\frac{{\partial {F_e}}}{{\partial H_{yi}^e}} = \sum\limits_{j = 1}^{ne} {\int\limits_{e}{{\left[ { - \frac{{\partial N_i^e}}{{\partial x}}\frac{{\partial N_j^e}}{\partial {y}}H_{xj}^e + \left( {\frac{{\partial N_i^e}}{{\partial x}}\frac{{\partial N_j^e}}{{\partial x}} + \frac{{\partial N_i^e}}{{\partial z}}\frac{{\partial N_j^e}}{{\partial z}} - {k^2}N_i^eN_j^e} \right)H_{yj}^e - \frac{{\partial N_i^e}}{{\partial z}}\frac{{\partial N_j^e}}{\partial {y}}H_{zj}^e} \right]}}} dxdydz = \sum\limits_{j = 1}^{ne} {\int\limits_{e}{{{N_i}{N_j}dxdydz{R_{yj}}}}} ,}\\ {\frac{{\partial {F_e}}}{{\partial H_{zi}^e}} = \sum\limits_{j = 1}^{ne} {\int\limits_{e}{{ - \frac{{\partial N_i^e}}{{\partial x}}\frac{{\partial N_j^e}}{{\partial z}}H_{xj}^e - \frac{{\partial N_i^e}}{\partial {y}}\frac{{\partial N_j^e}}{{\partial z}}H_{yj}^e + \left[ {\left( {\frac{{\partial N_i^e}}{{\partial x}}\frac{{\partial N_j^e}}{{\partial x}} + \frac{{\partial N_i^e}}{\partial {y}}\frac{{\partial N_j^e}}{\partial {y}} - {k^2}N_i^eN_j^e} \right)H_{zj}^e} \right]}}} dxdydz = \sum\limits_{j = 1}^{ne} {\int\limits_{e}{{{N_i}{N_j}dxdydz{R_{zj}}}}},} \end{array}} \right.\end{equation}$$ (8) where |$N_i^e$| represents the shape function of the node i in the unit and ne indicates the number of nodes in the unit. Each unit sub-matrix will then be calculated to synthesise the overall equation, and the magnetic field value at each node can be obtained by solving the equation. 2.5. Equation solving and algorithm verifications The locations of the non-zero elements in the stiffness matrix, as well as the required storage space, were unknown. Therefore, if a statistical data structure had been used in the overall integration process, the storage space could not be reasonably allocated according to the characteristics of the stiffness matrix. Moreover, to ensure that all of the non-zero elements could be stored, it was necessary to allocate a large static space, which would result in a waste of storage space or even insufficient storage space. Therefore, to avoid the above-mentioned problems, a dynamic data structure-a double linked list was adopted in the present study. Following the assemblage of the total rigid matrix, the next problem to be faced was how to solve the large sparse linear complex-coefficient electromagnetic field equations that had been formed. The methods that are used to directly solve linear equations will determine the accuracy and speed of forward modelling using finite element methods. In the current study, the large-scaled sparse matrix parallel solver PARDISO was used for the solution process. In this study, after obtaining H, the magnetic induction component was calculated according to equation (9) as follows $$ \begin{equation}{B_T} = \mu H.\end{equation}$$ (9) Then, to verify the correctness of the finite element solution, the magnetic field generated by a loop source with l = 100 m in the uniform half space model with ρ = 100 Ω m−1 was calculated, and different frequency solutions of the center point of the loop were compared with the analytical solution. As detailed in figure 2, the finite element calculation results were basically consistent with the analytical solutions. Figure 2. Open in new tabDownload slide Comparison between the finite element solutions and analytical solutions of the uniform half-space model. (a) Real part and (b) imaginary part. Figure 2. Open in new tabDownload slide Comparison between the finite element solutions and analytical solutions of the uniform half-space model. (a) Real part and (b) imaginary part. 3. Ground NMR response calculations 3.1. Forward modelling of NMR Geomagnetic fields contain both the geomagnetic dip I and the geomagnetic declination D, and the vertical components of induced magnetic fields B⊥ are more complicated. In the present study, for the purpose of simplifying the calculations, the coordinate system was rotated. First, the axis z was the rotation axis, and the axes x and y were rotated clockwise by the angle D, which was then recorded as the coordinate system xʹyʹzʹ. Second, axis yʹ was the rotation axis, and the axes xʹ and zʹ were rotated clockwise by an angle I, which was then recorded as the coordinate system xʺyʺzʺ. Meanwhile, the relationship between the induced magnetic field |${B^{\prime\prime}_T}$| in the new coordinate system and the original coordinate system B⊥ was expressed as follows: $$\begin{equation}{B^{\prime\prime}_T} = {R_I}{R_D}{B_T},\end{equation}$$ (10) where $$\begin{equation*} {R_D} = \left[ {\begin{array}{@{}*{3}{c}@{}} {\cos D\,\,}&{\sin D\,\,}&0\\ { - \sin D\,\,\,\,\,}&{\cos D\,\,\,\,}&0\\ 0&0&1 \end{array}} \right],\quad {R_I} = \left[ {\begin{array}{@{}*{3}{c}@{}} {\cos I\,\,}&0&\,\,{\sin I}\\ 0&\,\,1&\,\,0\\ { - \sin I\,\,}&0&\,\,{\cos I} \end{array}} \right].\end{equation*}$$ Then, since the axis xʺ direction of |${B^{\prime\prime}_T}$| had coincided with the direction of the geomagnetic field B0, the vertical component was synthesised by the components of axis yʺ and axis zʺ as follows: $$\begin{equation}{B_ \bot } = \sqrt {{{({B^{\prime\prime}_y})}^2} + {{({{B^{\prime\prime}_z}})}^2}} .\end{equation}$$ (11) Similarly, when the normal skew αn and normal inclination βn of the coil changed, the induced magnetic field Bʹʺ after the coordinate rotation could be derived by rotating the coordinate system as follows: $$ \begin{equation} B^{\prime\prime\prime} = {R_\beta }{R_\alpha }{B_T}, \end{equation}$$ (12) where $$\begin{equation*} {R_\alpha } = \left[ {\begin{array}{@{}*{3}{c}@{}} {\cos {\alpha _n\,\,}}&{\sin {\alpha _n\,\,}}&0\\ { - \sin {\alpha _n\,\,}}&{\cos {\alpha _n\,\,}}&0\\ 0&0&1 \end{array}} \right],{R_\beta } = \left[ {\begin{array}{@{}*{3}{c}@{}} {\cos {\beta _n\,\,}}&0&\,\,{\sin {\beta _n}}\\ 0&1&0\\ { - \sin {\beta _n\,\,}}&0&\,\,{\cos {\beta _n}} \end{array}} \right].\end{equation*}$$ In the equation, Rα and Rβ are positive definite. Then, $$\begin{equation}{B_T} = R_\beta ^{ - 1}R_\alpha ^{ - 1}B^{\prime\prime\prime} = R_\beta ^{\rm{T}}R_\alpha ^{\rm{T}}B^{\prime\prime\prime}.\end{equation}$$ (13) The vector value Bʹʺ of the induced magnetic field in the coil coordinate system was then calculated. The four rotation matrixes could then be used to obtain |$B^{\prime\prime} = {R_I}{R_D}R_\beta ^TR_\alpha ^TB^{\prime\prime\prime}$| under the geomagnetic field coordinate system. Finally, the vertical component B⊥ of the induced magnetic field was calculated using equation (13). After obtaining the excitation field B⊥, the initial amplitude response of the NMR could be successfully calculated using equation (1), or the nuclear magnetic free decay signal of the given pulse moment could be calculated using equation (14) as follows: $$\begin{equation}V(q) = \iiint {K(q,r) \cdot n(r)dV},\end{equation}$$ (14) where V(q) represents the nuclear magnetic FID attenuation signal; K(q, r) = |$\frac{{{\omega _0}{M_0}}}{I}{B_ \bot }\sin ( {\gamma \frac{{{B_ \bot }}}{I}q} ) \cdot {e^{2i\xi (r)}}$|is the kernel function; n(r) denotes the water content and q indicates the pulse moment. 3.2. Layered formation NMR In the present study, for the purpose of verifying the forward program, the given geomagnetic field strength was set as 48 000 nT; magnetic dip angle was 60°; and the magnetic declination was 0°. The loop side of the emission rectangle was 10 m in length. Then, to design a layered model (shown as figure 3), two layers of medium were known resistivity ρ1 and ρ2, respectively. The thickness of the first layer was h1, and the second layer of medium contained a 40 × 40 × 15 m rectangular water body. The water content was 10% and the water body resistivity was 100 Ω m−1. Then, since the known resistivity of the first layer ρ1 was 10 Ω m−1 and that of the second layer was ρ2 = 500 Ω m−1, the thicknesses h1 of the first layer were changed to 10, 20, 40 and 60 m, respectively. The initial amplitude response curve of the NMR is shown in figure 4. Figure 3. Open in new tabDownload slide Change the upper cover thickness to set a series of water-bearing bodies and calculate the initial amplitude response curve of the NMR. Figure 3. Open in new tabDownload slide Change the upper cover thickness to set a series of water-bearing bodies and calculate the initial amplitude response curve of the NMR. Figure 4. Open in new tabDownload slide NMR responses of subsurface water-bearing body with the first conductive layer of different thicknesses. Depth variation in the legend corresponds to different thicknesses of the first conductive layer. Figure 4. Open in new tabDownload slide NMR responses of subsurface water-bearing body with the first conductive layer of different thicknesses. Depth variation in the legend corresponds to different thicknesses of the first conductive layer. 4. Complicated emission source NMR responses 4.1. Eight-shaped loop NMR responses In this study, a given rectangular coil with a side length of 10 m was placed directly above the model. The model resistivity parameters were same as those of the layer model shown as figure 3. The water content was 10%, and buried depths were 10, 20, 40 and 60 m, respectively. The NMR responses are shown in figure 5a. The shape of the coil was changed to an ‘8’ shape to keep the emission magnetic moment unchanged, and its NMR responses are shown in figure 5b. Figure 5. Open in new tabDownload slide Initial amplitude curves of the different emission liner pulse moments. (a) Loop-shaped coil emission and (b) 8-shaped coil emission. Figure 5. Open in new tabDownload slide Initial amplitude curves of the different emission liner pulse moments. (a) Loop-shaped coil emission and (b) 8-shaped coil emission. It can be seen in figure 5 that the NMR responses excited by the 8-shaped loop under the same emission magnetic moment were stronger than those of the loop device. Therefore, under the same circumstances, the 8-shaped excitation mode displayed a certain anti-noise property. 4.2. Three-dimensional NMR simulations of complex field sources in undulating terrain Using the method described in this research study, the NMR imaging experiments involving the surfaces of karst water in a mountainous area were simulated to examine the effects of undulating topography, complex shape emission sources, and non-uniform resistivity. The model detailed in figure 6 shows a karst channel and a karst water storage model, and the surface is a real undulating mountain terrain. The underground areas contained a water storage cave and two confluent karst channels. The surface area was covered by Quaternary loess with a resistivity of 100 Ω m−1, and the lower part of loess layer was the limestone with a resistivity of 500 Ω m−1. The resistivity of the karst channel and the water-filled cave were both 30 Ω m−1. Figure 6 shows the fluctuations on the surfaces and the coil placement and illustrates the discrete case of the numerical calculations. The magnetic field distributions of the surface-emitting coil were calculated using an algorithm, and the magnetic field was distorted due to the presence of a surface low-resistance body. Figure 6. Open in new tabDownload slide Karst channel and a karst water storage model with real mountain terrain. Figure 6. Open in new tabDownload slide Karst channel and a karst water storage model with real mountain terrain. After obtaining the underground magnetic field, the kernel function K(q, r) was calculated. The kernel function distribution is shown in figure 7a and b. Compared with the uniform half-space NMR integral kernel function, the kernel function in the graph is affected not only by the low resistance region, but also by the topographic fluctuation, especially in the area close to the surface and the return line source, the magnetic field is seriously affected by topography, which leads to the complex change of the kernel function. The different pulse moment FID responses were calculated based on equation (13), and the initial amplitude was extracted. The results are detailed in figure 8. Figure 7. Open in new tabDownload slide The distribution of kernel function K(q, r) in the mountainous area. (a) Three-dimensional distribution map and (b) XOZ profile slice. Figure 7. Open in new tabDownload slide The distribution of kernel function K(q, r) in the mountainous area. (a) Three-dimensional distribution map and (b) XOZ profile slice. Figure 8. Open in new tabDownload slide Initial amplitude and FID responses of the given model. Figure 8. Open in new tabDownload slide Initial amplitude and FID responses of the given model. 5. Conclusion The algorithm used in this study involved the loading of a thin wire source and the calculation of the total magnetic field using a scalar finite element L, which improved the loading method of the source. To apply the calculations of the excitation field of the different transmitting coils, the loading of the source was realised by adding a power current density term at the right end of the equation. The double rotation equation was then transformed into a scalar finite element equation, and the zero divergence condition was found to have not been satisfied. Therefore, a divergence condition was required to be imposed to avoid the problem of pseudo solutions. Then, after obtaining the excitation field distributions, the vertical components of the geomagnetic fields were obtained using matrix rotation. As a result, the three-dimensional forward modelling of the terrestrial NMR was eventually realised. It was found in this study that the positions and directions of the transmitting coils had major influences on the initial amplitudes. Therefore, terrain effects should be considered in areas with severe terrain fluctuations. It was observed that under the same magnetic moment conditions, the nuclear magnetic signals excited by the 8-shaped coil were slightly stronger than those of the loop. This conclusion had also been inferred from the distribution of the excitation field. Acknowledgements The authors would like to appreciably acknowledge the support for this study that was provided by the Natural Science Foundation of China (grant nos 2013CB03600, 41304114 and 51139004), and China postdoctoral science foundation projects (no. 2017M620475). Conflict of interest statement. None declared. References Avdeev D. , Knizhnik S. , 2009 . 3D integral equation modeling with a linear dependence on dimensions , Geophysics , 74 , F89 – F94 . Google Scholar Crossref Search ADS WorldCat Braun M. , Hertrich M. , Yaramanci U. , 2005 . Study on complex inversion of magnetic resonance sounding signals , Near Surface Geophysics , 3 , 155 – 163 . Google Scholar Crossref Search ADS WorldCat Braun M. , Yaramanci U. , 2008 . Inversion of resistivity in magnetic resonance sounding , Journal of Applied Geophysics , 66 , 151 – 164 . Google Scholar Crossref Search ADS WorldCat Cai H. , Hu X. , Li J. , 2017a . Parallelized 3D CSEM modeling using edge-based finite element with total field formulation and unstructured mesh , Computers & Geosciences , 99 , 125 – 134 . Google Scholar Crossref Search ADS WorldCat Cai H. , Hu X. , Xiong B. , 2017b . Finite element time domain modeling of controlled-source electromagnetic data with a hybrid boundary condition , Journal of Applied Geophysics , 145 , 133 – 143 . Google Scholar Crossref Search ADS WorldCat Cui T. , Chew W. , 2000 . Modeling of arbitrary wire antennas above ground , IEEE Transactions on Geoscience and Remote Sensing , 38 , 357 – 365 . Google Scholar Crossref Search ADS WorldCat Haber E. , 1999 . Modeling 3D EM using potentials and mixed finite elements, three dimensional electromagnetics , SEG Geophysical Developments Series , 7 , 12 – 15 . WorldCat Haber E. , Ascher U.M. , 2001 . Fast finite volume simulation of 3D electromagnetic problems with highly discontinuous coefficients , SIAM Journal on Scientific Computing , 22 , 1943 – 1961 . Google Scholar Crossref Search ADS WorldCat Haber E. , Heldmann S. , 2007 . An octree multigrid method for quasi-static Maxwell's equations with highly discontinuous coefficients , Journal of Computational Physics , 223 , 783 – 796 . Google Scholar Crossref Search ADS WorldCat Haber E. , Heldmann S. , Ascher U. , 2007 . Adaptive finite volume method for distributed non-smooth parameter identification , Inverse Problems , 23 , 1659 . Google Scholar Crossref Search ADS WorldCat Hertrich M. , Braun M. , Yaramanci U. , 2005 . Magnetic resonance soundings with separated transmitter and receiver loops , Near Surface Geophysics , 3 , 141 – 154 . Google Scholar Crossref Search ADS WorldCat Hertrich M. , Yaramanci U. , 2002 . Joint inversion of surface nuclear magnetic resonance and vertical electrical sounding , Journal of Applied Geophysics , 50 , 179 – 191 . Google Scholar Crossref Search ADS WorldCat Jahandari H. , Farquharson C.G. , 2014 . A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids , Geophysics , 79 , E287 – E302 . Google Scholar Crossref Search ADS WorldCat Jaysaval P. , Shantsev D. , de la Kethulle de Ryhove S. , 2014 . Fast multimodel finite-difference controlled-source electromagnetic simulations based on a Schur complement approach , Geophysics , 79 , E315 – E327 . Google Scholar Crossref Search ADS WorldCat Lehmann-Horn J.A. , Hertrich M. , Greenhalgh S.A. , Green A.G. , 2011 . Three-dimensional magnetic field and NMR sensitivity computations incorporating conductivity anomalies and variable-surface topography , IEEE Transactions on Geoscience and Remote Sensing , 49 , 3878 – 3891 Google Scholar Crossref Search ADS WorldCat Legchenko A. , Ezersky M. , Camerlynck C. , 2009 . Joint use of TEM and MRS methods in a complex geological setting , Comptes Rendus Geoscience , 341 , 908 – 917 . Google Scholar Crossref Search ADS WorldCat Legchenko A. , Ezersky M. , Girard J. , 2008 . Interpretation of magnetic resonance soundings in rocks with high electrical conductivity , Journal of Applied Geophysics , 66 , 118 – 127 . Google Scholar Crossref Search ADS WorldCat Li G. , Han B. , 2017 . Application of the perfectly matched layer in 2.5D marine controlled-source electromagnetic modeling , Physics of the Earth and Planetary Interiors , 270 , 157 – 167 . Google Scholar Crossref Search ADS WorldCat Li J. , Hu X. , Zeng S. , 2013 . Three-dimensional forward calculation for loop source transient electromagnetic method based on electric field Helmholtz equation , Chinese Journal Of Geophysics , 56 , 4256 – 4267 . WorldCat Li J. , Lu X. , Farquharson C.G. , Hu X. , 2018 . A finite-element time-domain forward solver for electromagnetic methods with complex-shaped loop sources , Geophysics , 83 , E117 – E132 . Google Scholar Crossref Search ADS WorldCat Li J. , Zhu Z. , Liu S. , 2011 . 3D numerical simulation for the transient electromagnetic field excited by the central loop base on the vector finite-element method , Journal of Geophysics Engineering , 8 , 560 – 567 . Google Scholar Crossref Search ADS WorldCat Lin J. , 2010 . Situation and progress of nuclear magnetic resonance technique for groundwater investigations , Progress in Geophysics , 25 , 681 – 691 . WorldCat Mogi T. , 1996 . Three-dimensional modeling of magnetotelluric data using finite element method , Journal of Applied Geophysics , 35 , 185 – 189 . Google Scholar Crossref Search ADS WorldCat Ryu J. , Morrison H.F. , Ward S.H. , 1970 . Electromagnetic fields about a loop source of current , Geophysics , 35 , 862 – 896 . Google Scholar Crossref Search ADS WorldCat Newman G.A. , Alumbaugh D.L. , 1995 . Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences , Geophysical Prospecting , 43 , 1021 – 1042 . Google Scholar Crossref Search ADS WorldCat Raiche A.P. , 1987 . Transient electromagnetic field computations for polygonal loops on layered earths , Geophysics , 52 , 785 – 793 . Google Scholar Crossref Search ADS WorldCat Shushakov O. , 1996 . Groundwater NMR in conductive water , Geophysics , 61 , 998 – 1006 . Google Scholar Crossref Search ADS WorldCat Wan L. , Lin T. , Lin J. , 2013 . Joint inversion of MRS and TEM data based on adaptive genetic algorithm , Chinese Journal of Geophysics , 56 , 3728 – 3740 . WorldCat Wannamaker P.E. , Hohmann G.W. , SanFilipo W.A. , 1984 . Electromagnetic modeling of three-dimensional bodies in layered earths using integral equations , Geophysics , 49 , 60 – 74 . Google Scholar Crossref Search ADS WorldCat Wong A. , Wang X. , Liu G. , 2007 . Nonlinear inversion of surface nuclear magnetic resonance over electrically conductive medium , Chinese Journal of Geophysics , 50 , 890 – 896 . WorldCat Xiong Z. , 1992 . Electromagnetic modeling of 3-D structures by the method of system iteration using integral equations , Geophysics , 57 , 1556 – 1561 . Google Scholar Crossref Search ADS WorldCat Zhdanov M. , Dmitriev V.I. , Gribenko A. , 2007 . Integral electric current method in 3-D electromagnetic modeling for large conductivity contrast , IEEE Transactions on Geoscience and Remote Sensing , 45 , 1282 – 1290 . Google Scholar Crossref Search ADS WorldCat Zhou N. , Xue G. , Li M. , 2011 . Response of large Loop TEM based on approximation of electric dipole , Coal Geology and Exploration , 39 , 49 – 54 . WorldCat © The Author(s) 2019. Published by Oxford University Press on behalf of the Sinopec Geophysical Research Institute. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
TI - The three-dimensional finite element forward modelling of complex excitation source NMR
JF - Journal of Geophysics and Engineering
DO - 10.1093/jge/gxz096
DA - 2020-02-01
UR - https://www.deepdyve.com/lp/oxford-university-press/the-three-dimensional-finite-element-forward-modelling-of-complex-JrI8oPRup0
SP - 127
VL - 17
IS - 1
DP - DeepDyve
ER -