Analysis of Acoustic–Magnetic Fields Induced by Underwater Pressure Wave in a Finite-Depth Ocean

Abstract

As underwater disturbances (natural or artificial) occur in the ocean, moving seawater crossing the geomagnetic fields will produce weak circular currents. These currents can induce measurable magnetic fields, which might be useful for monitoring ocean internal waves using aeromagnetic survey. In this research, a spectral-element method (SEM) based on Gauss–Lobatto–Legendre (GLL) polynomials is presented to characterize the magnetic field induced by the underwater pressure waves. A concise mathematical model is established through combining the acoustic wave equations and Maxwell’s equations. Specifically, the acoustic–magnetic coupling simulation adopts the nodal-based SEM for acoustic analysis and edge-based SEM for electromagnetic analysis. The proposed SEM has spectral accuracy, as the error exponentially decreases with the order of the basis functions. Additionally, by adopting an independent modeling and mesh scheme in two solvers, respectively, the waste of computing resources is avoided. The experimental analysis demonstrates that the induced magnetic fields mechanically propagate with the acoustic wave, producing the pseudo-radiation phenomenon. The signals of these magnetic fields may extend for tens of kilometers and exist for hours under certain circumstances, which provide a theoretical basis for underwater target identification via high-sensitivity atomic magnetometer.

1. Introduction

Large-scale underwater disturbances are usually caused by the internal waves that are widely found in the deep ocean [1]. The nonlinear mechanisms of internal wave generation are complicated, mainly including tidal internal waves and source-induced internal waves. The former has an important influence on marine structures and is an environmental factor that must be considered in the design of marine engineering structures, while the latter caused by oscillations of underwater pressure waves may be used for purposes such as monitoring underwater targets [2]. Since seawater is a weak dielectric, moving water particles caused by underwater disturbances crossing the geomagnetic field will create weak induced currents, which in turn generate induced magnetic fields in and around the ocean [3,4]. Therefore, the establishment of an acoustic–magnetic coupling mathematical model may provide theoretical reference for aeromagnetic anomaly detection and underwater target identification [5].

During the past few decades, several approaches have been developed to evaluate the marine magnetic fields [6,7,8,9]. The standard approach for calculating marine magnetic signals caused by source-induced waves usually involves two essential steps. The first step is to establish the velocity model of seawater according to the dynamic mechanism, and the second step is to simulate the weak induced magnetic field based on the velocity model and Maxwell’s equations. The literature [10] proposed a numerical method for electromagnetic waves caused by a moving underwater target in layered seawater. The water velocity is first computed using the FVM (Finite volume method). Then, the electromagnetic wave is simulated using the FEM (Finite element method). Similarly, a study in [11] provided a common calculation method of electromagnetic waves induced by moving water in an infinite-depth ocean through the geomagnetic field. Based on this approach, two kinds of motions were investigated. Correspondingly, [12] analyzed the magnetic-field disturbance caused by the movement of a spherical dielectric in a conducting fluid. The solution for this problem is presented in the case that the magnetic field can be in any direction and the fluid is laminar around this ball. Recent research in [13] proposed a theoretical formula for calculating the induced magnetic field and studied the influence of environmental variables and geometric features of underwater moving targets on their behavior. Conclusively, previous studies [9,10,11] of induced magnetic field radiation phenomena demonstrated that the magnitudes of induced fields at the sea surface attain several nano-Tesla (nT). Therefore, the optical-pumping atomic magnetometers [14,15] with pico-Tesla (pT) precision may be utilized to accurately distinguish the signal from submerged moving target out of ambient noise.

For the investigation of a marine magnetic field, the difficulty of theoretical research lies in the establishment of accurate flow velocity models corresponding to different motion states of seawater. Since we focus on the source-induced internal waves here, the difficulty may be resolved by describing the velocity model as induced by the superposition of submerged acoustic sources. Therefore, this research investigates the multi-physics simulation technique [16,17,18] based on the spectral element method (SEM) for an aeromagnetic survey, as shown in Figure 1. Specifically, we first employ the nodal-based SEM to solve the Helmholtz equation for the acoustic pressure to obtain the flow velocity distribution in the ocean. Since the acoustic wave equations and Maxwell’s equations are coupled by flow velocity, we can further analyze the vector Helmholtz equation for a magnetic field by introducing the edge-based SEM, which has spectral accuracy and can solve large-scale problems with relatively few unknowns [19,20]. Utilizing the GLL (Gauss–Lobatto–Legendre) points for designing the quadrature integration points, the SEM exhibits a better convergence rate over traditional FEM [21,22,23].

The contribution of this thesis can be summarized as: (1) A concise mathematical model of underwater acoustic–magnetic coupling phenomenon is established, which may be an alternative solution to the existing theoretical difficulties. (2) The hybrid SEMs with scalar and vector basis functions were developed to construct independent modeling, which have better calculation efficiency than traditional FEMs for solving large-scale simulation problems in an aeromagnetic survey. (3) The acoustic–magnetic coupling phenomenon is analyzed, the magnitudes of induced magnetic fields above the sea are evaluated, and visual representations of the pseudo-radiation field are reported.

2. Acoustic–Magnetic Coupling Formula

To characterize the acoustic–magnetic coupling phenomenon [24,25,26] induced by underwater pressure waves, a multi-physics model can be established through combining the acoustic wave equations and Maxwell’s equations. For deriving the acoustic wave equation of interest, two fundamental laws (balance of momentum and conservation of mass) are employed [27,28]. Commonly, problems of linear acoustics refer to small perturbations of ambient quantities. The small fluctuating parts of pressure, density, and flow velocity are represented as $\tilde{p}$, $\tilde{\rho }$, and $\tilde{v}$, respectively. For the sake of simplicity, we do not take ambient flows into account, which allows for the governing equation of the acoustic wave to be written as

$$\left[\begin{array}{ccc}\partial t& {\rho }_0▽·& 0\\ 0& {\rho }_0\partial t& ▽\\ -c_s^2& 0& 1\end{array}\right]\left[\begin{array}{c}\tilde{\rho }\\ \tilde{v}\\ \tilde{p}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]$$

(1)

where $c_s$ is the speed of sound in fluids, and ${\rho }_0$ is the ambient quantity of the density. Reducing the above problem to a single variable helps further description. The variable of interest is the pressure fluctuation, hereafter referred to as the acoustic pressure. The local conservation of mass, the Euler equation as the balance of momentum, and the constitutive equation are all summarized into one partial differential equation, i.e., the wave equation. For time–harmonic problems, we have the time-dependence $\tilde{p}(x,t)=p\left(x\right)e^{-j\omega t}$; the acoustic-pressure Helmholtz equation can be described as

$${▽}^{2}p\left(x\right)+k^{2}p\left(x\right)=0x\in \mathrm{\Omega }\subset {\mathbb{R}}^3.$$

(2)

Here, $j=\sqrt{-1}$ denotes the imaginary unit, $k=\omega /c_s$ represents the wave number with the circular frequency $\omega$, $\mathrm{\Omega }$ is the computational domain with closed boundary $\partial \mathrm{\Omega }$. To further obtain a weak-form equation, we test (2) with the testing function $\tau \left(x\right)$ to provide the weak form

$${\int }_{\mathrm{\Omega }}[▽\tau \left(x\right)·▽p\left(x\right)-k^2\tau \left(x\right)p\left(x\right)]dv-{\int }_{\partial \mathrm{\Omega }}\tau \left(x\right)\frac{\partial p\left(x\right)}{\partial n}ds=0$$

(3)

where n represents the outward normal vector on the boundary surface and $\partial /\partial n$ is the normal derivative.

To increase the higher order accuracy and simulation efficiency, we employ a nodal-based GLL basis function to solve the above equation. The scalar basis function can be written in the 3-D reference domain $(\xi ,\eta ,\zeta )$ as

$${\widehat{\mathrm{\Psi }}}_{rst}={\varphi }_r^N\left(\xi \right){\varphi }_s^N\left(\eta \right){\varphi }_t^N\left(\zeta \right)$$

(4)

Here, ${\varphi }_i^N$ is the Nth-order GLL basis functions in the 1-D reference domain [29], N is the interpolation order, and i = 0, 1, …, N. Consequently, the discrete approximation for $p\left(x\right)$ in reference domain can be expressed as

$$\widehat{p}(\xi ,\eta ,\zeta )=\sum _{r=0}^N\sum _{s=0}^N\sum _{t=0}^N{\widehat{k}}_{rst}{\widehat{\mathrm{\Psi }}}_{rst}$$

(5)

where ${\widehat{k}}_{rst}$ represents the degrees of freedom (DOF) at the nodal point $({\xi }_r,{\eta }_s,{\zeta }_t)$. Finally, the covariant mappings from the reference domain to the physical domain can be expressed as

$$\mathrm{\Psi }=\widehat{\mathrm{\Psi }} \mathrm{and} ▽\mathrm{\Psi }={\left(JM\right)}^{-1}\widehat{▽}\widehat{\mathrm{\Psi }}$$

where $\mathrm{\Psi }$ denotes the scalar basis function in the physical domain, and JM represents the Jacobian matrix [30]. On top of that, appropriate boundary conditions were imposed to describe the horizontally stratified ocean, as will be discussed later. Once the pressure $p\left(x\right)$ is obtained, the flow velocity vector $\mathit{v}\left(x\right)$ in the fluid can be calculated according to

$$\mathit{v}\left(x\right)=▽p\left(x\right)/\left(j\omega \rho \right)$$

(6)

After the calculation of the flow velocity distribution in the ocean, the acoustic and magnetic problems could be coupled via Ohm’s law for fields in motion.

On the other hand, the time-domain Maxwell’s equations in inhomogeneous media with permittivity $ϵ$, permeability $\mu$, and conductivity $\sigma$ can be expressed as

$$\left[\begin{array}{cc}-▽\times & -\mu \partial t\\ -ϵ\partial t-\sigma & ▽\times \end{array}\right]\left[\begin{array}{c}\mathit{E}\\ \mathit{H}\end{array}\right]=\left[\begin{array}{c}{\mathit{M}}_{\mathit{i}}\\ {\mathit{J}}_{\mathit{i}}\end{array}\right]$$

(7)

where ${\mathit{M}}_{\mathit{i}}$ and ${\mathit{J}}_{\mathit{i}}$ are the impressed magnetic and electric current densities, respectively.

In magneto-hydrodynamics (MHD), we consider only the magnitude of velocity $\left|\mathit{v}\right|\ll c_l$ ($c_l$ represents the speed of light), which means that the contribution of the time derivative of $ϵ\mathit{E}$ is much smaller than the current density. Therefore, we can omit the dependence of the Ampere–Maxwell law on $\partial \mathit{D}/\partial t$ and use the pre-Maxwell form of $▽\times \mathit{H}={\mathit{J}}_{\mathbf{c}}+{\mathit{J}}_{\mathit{i}}$. Furthermore, since the water particles are moving along a velocity field $\mathit{v}$, the conducted current density ${\mathit{J}}_{\mathbf{c}}$ is additionally related to the magnetic flux density $\mathit{B}$, and the Ohm’s law for the fields in motion becomes ${\mathit{J}}_{\mathbf{c}}=\sigma (\mathit{E}+\mathit{v}\times \mathit{B})$. Note that the impressed current density ${\mathit{J}}_{\mathit{i}}=0$ in this case; then, the equation describing the evolution of $\mathit{B}$ in an electrically conducting fluid can be expressed as

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}-▽\times (\mathit{v}\times \mathit{B})+\lambda ▽\times (▽\times \mathit{B})=0\end{array}$$

(8a)

$$\begin{array}{c}\hfill ▽·\mathit{B}=0\end{array}$$

(8b)

where $\lambda ={\left(\sigma \mu \right)}^{-1}$ is the magnetic diffusivity coefficient. If we further express the current density ${\mathit{J}}_{\mathbf{s}}=\sigma (\mathit{v}\times {\mathit{B}}_E)$, which in this case is the electric field induced by the movement of the water particles crossing the geomagnetic field ${\mathit{B}}_E$, we arrive at

$${▽}^2\mathit{A}-j\omega {\lambda }^{-1}\mathit{A}=-\mu {\mathit{J}}_{\mathbf{s}}$$

(9)

Here, $\mathit{A}$ represents the magnetic vector potential defined through $▽\times \mathit{A}=\mathit{B}$. By testing (9) with function $\mathit{w}\left(x\right)$, we obtain another weak-form equation related to magnetic fields:

$$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}[▽\mathit{w}\left(x\right)·▽\mathit{A}\left(x\right)+j\omega {\lambda }^{-1}\mathit{w}\left(x\right)\mathit{A}\left(x\right)]dv\\ \hfill -{\int }_{\partial \mathrm{\Omega }}\mathit{w}\left(x\right)\frac{\partial \mathit{A}\left(x\right)}{\partial n}ds=\mu {\int }_{\mathrm{\Omega }}\mathit{w}\left(x\right){\mathit{J}}_{\mathbf{s}}\left(x\right)dv\\ \hfill \end{array}$$

(10)

For solving the vector Helmholtz equation defined by (10), the spectral element discretization is employed, and the mixed-order curl-conforming vector basis functions are introduced as follows:

$${\widehat{\mathrm{\Phi }}}_{rst}^{\xi }=\widehat{\xi }{\varphi }_r^{N-1}\left(\xi \right){\varphi }_s^N\left(\eta \right){\varphi }_t^N\left(\zeta \right){\widehat{\mathrm{\Phi }}}_{rst}^{\eta }=\widehat{\eta }{\varphi }_r^N\left(\xi \right){\varphi }_s^{N-1}\left(\eta \right){\varphi }_t^N\left(\zeta \right){\widehat{\mathrm{\Phi }}}_{rst}^{\zeta }=\widehat{\zeta }{\varphi }_r^N\left(\xi \right){\varphi }_s^N\left(\eta \right){\varphi }_t^{N-1}\left(\zeta \right)$$

which are parallel to the $\xi -$, $\eta -$, and $\zeta -$ axes in the 3-D reference domain, respectively. Thus, the discrete approximation for $\mathit{A}$ in the reference domain $\widehat{\mathit{A}}=(A_{\xi },A_{\eta },A_{\zeta })$ is provided by

$$\begin{array}{c}\hfill \widehat{\mathit{A}}(\xi ,\eta ,\zeta )=\sum _{r=0}^{N-1}\sum _{s=0}^N\sum _{t=0}^N{\widehat{h}}_{rst}^{\xi }{\widehat{\mathrm{\Phi }}}_{rst}^{\xi }+\sum _{r=0}^N\sum _{s=0}^{N-1}\sum _{t=0}^N{\widehat{h}}_{rst}^{\eta }{\widehat{\mathrm{\Phi }}}_{rst}^{\eta }+\sum _{r=0}^N\sum _{s=0}^N\sum _{t=0}^{N-1}{\widehat{h}}_{rst}^{\zeta }{\widehat{\mathrm{\Phi }}}_{rst}^{\zeta }\end{array}$$

(11)

where ${\widehat{h}}_{rst}^{\xi }$, ${\widehat{h}}_{rst}^{\eta }$, and ${\widehat{h}}_{rst}^{\zeta }$ denote the DOFs at the nodal point $({\xi }_r,{\eta }_s,{\zeta }_t)$. Finally, the covariant mappings from the reference domain to the physical domain can be expressed as

$$\mathrm{\Phi }={\left(JM\right)}^{-1}\widehat{\mathrm{\Phi }} \mathrm{and} ▽\times \mathrm{\Phi }=\frac{1}{\left|JM\right|}{JM}^T\widehat{▽}\times \widehat{\mathrm{\Phi }}$$

where $\mathrm{\Phi }$ denotes the vector basis function in the physical domain. Once the magnetic vector potential $\mathit{A}\left(x\right)$ is calculated, current rings may be defined and the magnetic flux density $\mathit{B}$ would provide a visual representation of the field distributions induced by the underwater pressure wave.

3. Numerical Experiment

For solving the aeromagnetic survey model, the multi-physics problem was split into acoustic and magnetic parts. For acoustic simulation in a finite-depth ocean (1000 m), the computational domain is shown in Figure 2. An acoustic monopole source with a frequency of 0.75 Hz inducing a pressure wave of 1 Pa is located at 100 m underwater. In order to accurately simulate the marine environment, free surface and rigid bottom boundary conditions are adopted to simulate the ocean surface and seabed, respectively. The perfectly matched layer (PML) is enforced on the rest of the computational domain to simulate the open boundary. The speed of sound in fluids is assumed to be $c_s=1500$ m/s.

The acoustic pressure and flow velocity distribution are calculated using the nodal-based SEM. Figure 2a demonstrates the pressure wave distribution within a range of 3 km. It can be seen that the source is located at the far left of the computational domain where the pressure is the strongest. For the far field, the pressure decays slowly as the distance increases. In addition, Figure 2b shows the radial component of the velocity distribution. We can see that the flow velocity is the largest near the source, and it decays slowly with the increasing distance. However, due to different boundary conditions on the upper and lower surfaces, the distribution of velocities on these two surfaces are obviously different; the same applies for pressure fields.

Figure 3a shows pressure values on a series of detectors along the radial direction (x-axis) at a depth of 500 m underwater. We can see that the pressure value is about ${10}^{-7}$ Pa at 50 km away from the source, which shows that pressure waves slowly decay in the ocean and can propagate a long distance. Furthermore, we compared the velocity values obtained by the nodal-based SEM with the commercial software COMSOL, as shown in Figure 3b. It can be seen that our result agrees well with the reference. The flow velocity slowly decays and the value is on the order of ${10}^{-11}$ m/s at 50 km away from the source, indicating that the velocity field of the pressure wave has a wide distribution range underwater, which will induce a wide range of magnetic fields, as will be discussed later.

For magnetic-field simulation in a finite-depth ocean, the computational domain $\mathrm{\Omega }{}^{\prime }$ surrounded by PML has a dimension of 50 km × 10 km. The grid size is 100 m and the total number of elements is about 50,000. As we employed 3-order basis functions in our simulation, this grid scale ensures computational accuracy. For the time iteration part, we adopted 0.005 s as the time step, and the total number of the time steps is 10,000. The other parameters are as follows: air layer $({ϵ}_{r,A}={\mu }_{r,A}=1)$ with a 6-km thickness, coastal seawater $({ϵ}_{r,SW}=81,{\mu }_{r,SW}=1$ and ${\sigma }_{SW}=4$ S/m) with a 1-km thickness, and the lithosphere $({ϵ}_{r,Lit}=12,{\mu }_{r,Lit}=1$ and ${\sigma }_{Lit}=0.005$ S/m) with a 3-km thickness. Since the sea level serves as the starting point for coordinates, the acoustic monopole source is located at (0, 0, and −100) m, as shown in Figure 4. The geomagnetic field ${\mathit{B}}_E=5\times {10}^4\left(nT\right)$ is assumed to be in the vertical direction.

The distribution of the induced magnetic field generated by the underwater pressure wave are shown in Figure 4a. It can be seen that due to the high conductivity of seawater, the induced magnetic field is stronger in the ocean region. Although the magnetic field decays with distance, there is no significant change at the far end. Moreover, Figure 4b demonstrates the vertical component of the magnetic field, indicating that the strongest field is distributed near the source. Note that, due to the huge difference in the conductivity between the ocean surface and seabed, a natural waveguide is formed in the ocean. The induced magnetic field is confined in the waveguide, in which it propagates, oscillates, and decays.

The edge-based SEM results of the induced magnetic fields along the receiver are compared with those obtained via COMSOL, as shown in Figure 5 and

Table 1. Calculation statistics of edge-based SEM and COMSOL.

Method/Software	DOFs	CPU Time (s)	Memory Consumption (MB)	Acceleration Technique
Edge-based SEM (3 Order)	176,525	716	6034	None
COMSOL (3 Order)	2,274,302	1605	36,712	Parallel (6 Core)

. Figure 5a shows the vertical component of the magnetic field on a series of receivers along the radial direction (x-axis) at 500 m above the sea level. We can see that it decays slowly in the air and propagates a long distance. On top of that, we compared the magnetic field amplitudes calculated using the edge-based SEM and COMSOL, as shown in Figure 5b. It can be seen that our result agrees well with the reference. The induced magnetic field decays rapidly in the first 5 km, but slowly in the later part.

Furthermore, the calculation statistics of the whole simulation process using SEM and COMSOL are summarized in Table 1. For the proposed method that employed 3-order basis functions, the degrees of freedom (DOF) is 176,525, memory consumption is 6034 MB, and the CPU time is 716 s. For the software COMSOL, which adopted 3-order FEM, the DOF is 2,274,302, memory consumption is 36,712 MB, and the CPU time is 1605 s. It is worth noting that COMSOL used 6-core parallel acceleration technology in the above simulation. The proposed technique is more efficient for solving the multi-physics problem due to the use of the SEM that has spectral accuracy with a small number of unknowns to analyze 3-D problems.

Finally, we discuss the flow velocity distributions (detectors located 300 m underwater) generated by pressure waves with different frequencies, different intensities, and different depths, and the induced magnetic fields on the magnetometers located 100 m above the sea surface. Figure 6a,b show the radial component of the velocity field and amplitude of magnetic field generated by the pressure wave of 1 Pa with a frequency of 0.5 Hz at a depth of 100 m underwater, respectively. It can be seen that the frequency of the velocity field is in good agreement with that of the pressure wave. However, due to the low intensity of the pressure wave, the amplitude of the induced magnetic field above the ocean is relatively small. Figure 6c,d demonstrate the radial component of the velocity field and amplitude of the magnetic field generated by the pressure wave of 1 kPa with a frequency of 1 Hz at a depth of 200 m underwater, respectively. We can see that the frequency of the velocity field is approximately 1Hz, which is in good agreement with frequency of the pressure wave. In addition, due to the high intensity of the pressure wave, the amplitude of the induced magnetic field over the sea surface is relatively large. Figure 6e,f show the radial component of the velocity field and that the amplitude of the magnetic field generated by the pressure wave of 1 MPa has a frequency of 10 Hz at a depth of 500 m underwater, respectively. It can be seen that, due to the high intensity of the pressure wave, the intensities of the velocity field and the induced magnetic field are significantly increased. It is worth noting that the wave curves become dense because of the high-frequency acoustic source. In general, as the frequency and intensity of the pressure wave increase, so do the frequency and intensity of the induced magnetic field, appearing as nonlinear characteristics.

4. Conclusions

This research discusses the mechanism of an induced magnetic field generated by underwater pressure waves and proposes a concise mathematical model to demonstrate the acoustic–magnetic coupling phenomenon. Since rapid simulations of large-scale electromagnetic exploration are still a challenge due to the massive computational burden, the simulation algorithms are carefully designed. Not only the nodal-based SEM is introduced to solve the acoustic pressure field but also the edge-based SEM is employed to construct the magnetic field. The accuracy of the proposed algorithm is verified by using the commercial software. The numerical results are provided to characterize the marine magnetic field induced by various pressure waves. It shows that the induction signals may extend for tens of kilometers, causing observations to be possible when using commercial magnetometers.