Remote Sensing of the Ocean Environment Using Finite Element Methods

Abstract

Oceans are a vast, complex world where underwater sound is the most efficient tool available to understand its detailed characteristics. However the underwater channel has a very complex geometrical and material structure and hence special techniques are required to model it. Analytical solutions are feasible only when one makes gross assumptions and approximations. Several numerical and semi-numerical techniques have been developed for estimating the sound field in the ocean channel. But no single method is capable of handling all possible environmental conditions, frequency, and ranges of interest in remote sensing problems. We explore in this chapter the scope and feasibility of finite element method in underwater remote sensing. The current study is based on a channel model with cylindrical symmetry and a time-harmonic source signal. A variational formulation is used to derive the finite element model for acoustical radiation, scattering and propagation in the ocean. A Bayliss-type radiation boundary condition is used to model the far field behaviour without the need to deal with a large solution domain. Since the ocean geometry can support several propagating, evanescent, and radiation modes, a penalty function approach is employed to impose the far field radiation condition. A distinct feature of the ocean channel is its depth-dependent sound speed. The eigensolution for this channel is required for imposing the radiation condition at the truncation boundary. We have cast this eigenproblem in a variational form and employed a Rayleigh-Ritz method to obtain an approximate eigensolution. This approach has provided a good approximation of the depth eigenmodes in a compact semi-analytic form. We have employed our finite element algorithm to model several range- and depth-dependent ocean problems. Our numerical study has established that our finite element algorithm gives accurate results with reasonable effort. In particular, our finite element approach is most appropriate for shallow water problems where the interaction of wave modes with irregular ocean bottom is quite complex. The penalty function approach employed to implement the radiation boundary condition has been found to be robust over a wide range of penalty scale factors. We have also extended this work for the case of irregular elastic sea bed. We continue to explore and further develop our finite element approach by applying it to several other ocean acoustic problems encountered in the remote sensing of ocean environment.

1. Introduction

Oceans are a vast, complex, mostly dark, optically opaque but acoustically transparent world which is only thinly sampled by today’s limited science and technology. Underwater sound

There exists a vast body of literature in remote sensing of ocean using electromagnetic and optical sensors from satellites. Although such methods have definite advantages in several aspects, they have serious limitations for sensing deep underwater channels.

is used as the premier tool to determine the detailed characteristics of physical and biological bodies and processes in the ocean. The distributions within the sea of the physical variables affect the transmission of sound. The wide range of acoustic frequencies and wavelengths, together with the diverse oceanographic phenomena that occur over full spectra of space and time scales, thus give rise to a number of interesting effects and opportunities. Because of its great practical importance, especially to naval submarine operations, ocean-acoustics research [1-5] has been driven by applications more than other branches of ocean science.

Acoustic remote sensing in a generic sense refers to sending out acoustic signals and recording the scattered waves, which is hence processed to ascertain the nature of target/obstruction that was encountered by the transmitted signal. This remote sensing in general involves transmission, processing of received signals and some form of inversion. This chapter is exclusively dedicated to accurate modeling of propagation and scattering of acoustic signals in the ocean channel.

The amplitude and phase of sound field generated by an acoustic source in the ocean can be deduced, in principle, by solving either the wave equation or the Helmholtz equation in the case of a harmonic acoustic source [1]. However, this procedure is generally difficult to implement due to the complexity of the ocean-acoustic environment: the sound-speed profile is usually non-uniform in depth and/or range, giving rise to waveguide focusing and shadowing effects; the sea surface is rough and time-dependent; the ocean floor is typically a very complex, rough boundary which may be inclined to the horizontal; and the bottom may be an elastic medium, capable of supporting shear waves along the ocean-bottom boundary. To compound the problem, various ocean processes, including internal waves and small-scale turbulence, introduce small fluctuations in the sound speed, which are responsible for significant acoustic fluctuations over long transmission paths.

Analytical solutions of the governing differential equations in underwater acoustics are not always feasible and can only be obtained if the sound speed of the water column and physical boundaries can be described in simple mathematical terms. This is rarely the case in reality and so it is generally necessary to employ approximate models. A variety of numerical techniques have been developed for estimating sound fields in the ocean, but no single method is capable of handling all possible environmental conditions, frequencies, and transmission ranges of interest in the applications. Even the existing ocean-acoustic propagation models [6-7] with restricted scope often take several hours to run on a supercomputer.

Several different approaches for the solution of the sound field in the ocean have evolved over the past few decades: ray tracing [8], normal-mode techniques [9] and coupled-mode models [10], the parabolic-equation approximation [11] and fast field programs (FFP) [12]. In this chapter, we will discuss in detail the scope of the finite element method [13-15] in ocean remote sensing applications. In order to motivate our finite element approach and put it in proper context, we briefly summarize the analytical and computational tools that are currently in use in the ocean remote sensing literature especially to point out their main merits and shortcomings.

2. Background

Ray-based methods [8-9] involve following the paths of a set of rays as they leave the source and tracking them as they propagate through the medium. They can be used for range-dependent and range-independent problems, but are most commonly used for range-independent problems. They are most useful for short-range, high-frequency modeling. Straightforward ray theory suffers from following drawbacks: (i) Need to deal with situations involving caustics and singularities. (ii) At each incidence on surface or bottom, each ray has to be “told” at what angle to go off, and with what percentage of total reflection. (iii) Since problems are almost entirely numerical, each variation is nearly as hard as the first try, e.g., a new source depth or a greater range. The main shortcoming of the ray method is the inherent high-frequency approximation.

A class of propagation models exist which gives the full-wave solution for the field in a horizontally stratified medium. This type of a model is known as “fast field program”. This technique is basically a numerical implementation of the integral transform technique for horizontally stratified media [12, 9]. The field solution is in the form of a wavenumber integral which is evaluated by numerical quadrature. This approach is distinguished by its use of the fast Fourier transform (FFT) to calculate the integral. FFPs determine the field which satisfies the Helmholtz equation or similar equations which include shear wave effects. The Helmholtz equation for the stratified medium is a partial differential equation in two independent variables, range and depth, and hence in principle could be solved by the application of two integral transforms [12]. For certain specific sound-speed profiles having a particular analytical form, this can be achieved, yielding an exact solution for the field. For a general sound-speed profile, however, the transform over depth is intractable and an alternative technique must be sought. Nevertheless, a transform over range can be applied, and this is the starting point of the FFP argument [16]. In contrast to the ray solution, the FFP model yields a result which is essentially exact. Starting from the Helmholtz equation for a stratified medium, the only additional approximation is that of using the asymptotic approximation to the Bessel function. This approximation turns out to include negligible errors beyond a wavelength or so from the source.

As an alternative to “exact” numerical propagation models, with their heavy computational overhead, a number of methods have been developed whose starting point is a parabolic equation [11]. Such an equation which is an approximation for the elliptic Helmholtz equation is valid over a small range of angles, usually, but not necessarily, extending about the horizontal. Given their inherently approximate nature, the parabolic-equation (PE) models are distinguished by a lack of precision, the extent of which depends by and large on the problem under consideration. They have acquired popularity amongst the ocean-acoustics community because they give the field over the entire water column with no additional effort and they can handle range-dependent environments. PE methods are often said to be valid within a cone of angles extending +/− 20° (narrow angle) and +/− 40° (wide angle) about the horizontal. One of the shortcomings of the PE models is that, when these angles are exceeded, the output continues to look reasonable, showing no obvious indication of error [7]. Apart from the excessive inaccuracy of these results, the lack of consistency among the PE codes highlights the general difficulty of assessing their performance in any given environment. Although the PE is relatively easy to implement, there is a price to be paid: a) it is valid over only a limited range of angles, a consequence of the paraxial approximation, and b) it is a one-way solution, capable of handling only outgoing waves, since incoming radiation, represented by a Hankel function of the second kind of zero order, is neglected in the solution. Little can be done to remedy the backscatter limitation, but considerable effort has gone into extending the angular range of the forward-scatter regime [17]. The advantage of the parabolic equation over the original Helmholtz equation is that the PE can be solved by a straightforward marching in range which requires much less computational effort. From a numerical point of view, this range marching is typically implemented using either standard finite difference techniques or using a fast Fourier transform as in the so-called split-step method. There are other approaches to solve the parabolic wave equation in ocean waveguides. Lee et al. [18] employed the finite difference method whereas Huang [19] used a finite element method to solve the PE.

The sound field in a horizontally stratified ocean can be expressed as an infinite sum of uncoupled normal modes plus one or more branch line integrals [9, 1]. At large ranges from sources, the branch line integral component is negligible and the field is given accurately by the normal-mode sum, but in the vicinity of the source, within the cycle distance of each mode, the integrals are significant and should be taken into account. If the environment shows some range dependence, through either the sound-speed profile or the boundary conditions, the field is no longer separable and strictly (uncoupled) normal-mode theory does not apply. However, provided the range dependence is sufficiently slow, the adiabatic approximation is valid, i.e. there is essentially no transfer of energy between modes as they propagate the channel. If the range dependence is too fast for the adiabatic approximation to hold, mode coupling is significant, which requires the calculation of the coupling coefficients—a time consuming procedure. The normal-mode method is typically accurate for ranges greater than the first 10 water depths or so, a figure which depends on the number of modes that are included in the solution. In the near field, more modes should be computed to closely predict the fields accurately. The normal-mode models tend to be thought of as providing solutions to range-independent problems. Range-dependent solutions can be obtained using (a) adiabatic mode theory or (b) coupled-mode theory. The later approach involves more computational cost but can provide more accurate results.

When the range dependence is too strong for mode coupling to be neglected, a different approach than the usual normal-mode theory is required. A complete two-way (i.e. including backscattering) solution to this problem has been formulated in terms of stepwise coupled normal modes [10]. The medium is sub-divided into a large number of thin vertical segments, in each of which the acoustic parameters are held constant in the range direction but are allowed to vary in depth. Across the segment boundaries, the pressure and horizontal particle velocity are required to be continuous. In this method, the field is expressed as a sum of local modes representing both outgoing and incoming cylindrical waves. Again, the modal eigenvalue problem has to be solved, and in this case, the Galerkin method is used, whereby the solution is expanded in a set of basis modes, yielding a tractable eigenvalue matrix problem [20]. This involves rather, complex coupling integrals which have to be evaluated for all modes at all segment boundaries. This method is computationally demanding, but it is essentially exact and forms the basis of the model [10, 21-22]. When the coupling effects are neglected, the full coupled-mode expressions reduce to the adiabatic approximation.

Ray tracing, normal-mode techniques, and coupled-mode models are accurate but computationally intensive; the parabolic equation is an approximation to the wave equation that has been solved using explicit and implicit finite difference schemes; Green’s function solutions (fast field programs) are essentially models for which exact solutions are available that cannot account for sound-speed variation. If the variation of sound-speed profile is independent of range, the ocean is said to be horizontally stratified. Several of the numerical ocean-acoustic propagation models assume horizontal stratification. The advantage, from the point of view of the computation, is that the solution field separates into range and depth components, which simplifies the calculation of the field considerably. The speed of sound in the ocean shows only small departures from 1500 m/s, but nevertheless its effect on sound propagation on the ocean is profound. In the deep ocean, for example, the profile acts as an acoustic waveguide, supporting propagation to long ranges with little attenuation. However, for a general ocean environment, which has a range-dependent sound-speed profile, an ocean bed having irregular geometry, and turbulence in the water column, none of the existing methods described above work satisfactorily.

3. Finite element method

For general ocean environments, the finite element method (FEM) [13-15] is a good choice for the numerical modeling of ocean-acoustic wave propagation because it is exact within the limits of numerical accuracy and can accurately account for all the scattering processes. Although the literature on finite element technique on wave scattering and propagation is extensive, the number of available FEM models for ocean remote sensing is fairly small [23-24]. Part of the reason for this is the large computational cost involved. However, we feel that for shallow-water applications, the FEM is both feasible and appropriate. The very nature of the waves to radiate into the far field when unbounded requires the domain to be truncated with an artificial boundary, on which an approximate radiation boundary condition should be imposed [25-29]. In the present work, a variational approach is used to derive the finite element approximation for time-harmonic acoustic wave propagation in an axisymmetric, heterogeneous oceanic waveguide, and a BGT-type boundary damper [26] is used to model the effect of the far field. Since a waveguide in general supports multiple propagating/radiation modes in the far field, a penalty function approach has been employed to impose the modal radiation boundary condition in conjunction with the orthogonality property of the depth modes of the waveguide.

In our finite element model for depth- and range-dependent waveguides, the eigensolution of the depth problem is required for the imposition of the radiation condition at the truncation boundary. Unfortunately, the depth eigenproblem could be solved exactly only for a few special profiles. In view of this, several numerical methods have been developed to solve the depth problem [1]. Porter and Reiss [30] employ a finite difference model for the depth equation, and the resulting algebraic eigenproblem has been solved using a combination of iterative techniques and Richardson extrapolation to obtain the radial wavenumbers and modal vectors to a great degree of precision. For our finite element model [31-32], it would be convenient to have the depth modes in a compact analytical form. We have accomplished this by adopting the following procedure: The depth eigenproblem is cast in a variational form by suitably defining a functional. The classical Rayleigh–Ritz (RR) method is employed to find a variational approximation to the eigensolution of the depth problem in ocean-acoustic waveguides. The depth modes thus obtained have a compact semi-analytical form in contrast to methods using finite difference or other finite element methods. An interesting feature of the model is that the trial functions are derived from an isovelocity problem that has an exact solution. It is important to note that such trial functions automatically satisfy even the dynamic interface condition at the seabed, thus contributing to the accuracy of the numerical model. Our procedure has been tested for several different ocean profiles and the results compare well against those obtained using the method of Porter and Reis [33]. The proposed model thus provides an accurate representation of the depth eigenmodes in a compact semi-analytical form.

We have chosen several isovelocity waveguide examples, for which analytical solutions are available, to validate the FE model developed and ascertain its versatility to impose modal radiation boundary condition. We have confirmed the efficacy of the FE model by applying it to several examples of depth- and range-dependent waveguides. This numerical study establishes that our FE model gives accurate results with reasonable computational effort. The penalty function approach employed to implement the radiation boundary condition has been found to be robust over a wide range of penalty scale factors. We have also extended this work for the case of irregular elastic seabed. We continue to explore and further develop our FE model by applying it to several other ocean-acoustic problems encountered in the remote sensing of ocean environment.

4. Governing equations and boundary conditions

The fluid domain $$\mathrm{\Omega }={\mathrm{\Omega }}_I+{\mathrm{\Omega }}_O$$ (Fig. 1) of the waveguide problem consists of the inner domain $${\mathrm{\Omega }}_I$$ truncated by the artificial radiation boundary $$S_R$$, and the outer domain $${\mathrm{\Omega }}_O$$ (far-field domain). The waveguide is assumed to be axially symmetric about the vertical axis containing a source at depth $$z_s$$, with $$r$$ denoting the radial coordinate or the range. It is bounded at the top by the $$z=0$$ plane, which is the air–sea interface ($$S_F$$), and at the bottom by a seabed of arbitrary topography ($$S_B$$). The waveguide is assumed to have unbounded range. For time-harmonic linear acoustic waves with the pressure field denoted as $$\widehat{p}(r,z,t)=p(r,z)e^{-i\omega t}$$, $$\omega$$ being the circular frequency of the source, the governing equation is given by

$$\rho \nabla \left(\frac{1}{\rho }\nabla p\right)+k^{2}p=-\frac{1}{2\pi r}{f}_o\delta (r)\delta \left(z-z_s\right),$$E1

where $$\nabla$$ is the gradient operator, $$\rho$$ the density of the acoustic fluid, k the acoustic wavenumber, c the local speed of sound, and $$f_o$$ defines the point source at $$r=0$$ and $$z=z_s$$.

Figure 1.

Geometry of the ocean waveguide

Considering the large impedance mismatch between air and water, a pressure release boundary condition may be used at the free surface. Thus,

$$p=0\text{on}{S}_F$$E2

As the waves encounter the seabed, there is partial reflection and the remaining energy is transmitted into the seabed. A part of the transmitted waves may be coupled back into the water column because of refraction through the sediment layers. However, for now, a rigid bottom is assumed, for which the normal derivative of the pressure should vanish at the bottom boundary. In other words,

$$\frac{\partial p}{\partial n}=0\text{on}{S}_{\text{B}},$$E3

where $$S_{\text{B}}$$ denotes the sea bottom.

For the purpose of FE modeling, the waveguide, which is unbounded in range, is truncated at $$r=r_b$$, and the truncation boundary is treated as the radiation boundary $$S_R$$, on which a suitable approximate radiation condition should be imposed. Here, the boundary damper approach [26] has been adopted. The first-order cylindrical damper equation may be written as

$$\frac{\partial p_m}{\partial n}+\alpha p_m=0,m=1,2,\cdots ,M\text{on }{S}_R,$$E4

where M denotes the number of propagating modes, and the damper coefficient $${\alpha }_m$$ associated with the m-th mode is given by

$${\alpha }_m=\frac{1}{2r}-i{k}_{rm},$$E5

where $$k_{rm}$$ denotes a horizontal wavenumber. It may be noted that Eq. (5) is exact for the asymptotic form of a cylindrically symmetric wave. On the truncation boundary $$S_R$$, acoustic pressure may be expressed as a sum of normal pressure modes as

$$p(z)={\displaystyle \sum _{m=1}^M{p}_m(z)},$$E6

where $$p(z)$$, m=1,2,..., are the normal modes of propagation for the problem in Eq. (1). Following Fix and Marin [31], the radiation boundary condition for the waveguide problem may be written, using Eqs. (4) and (6), as

$$\frac{\partial p}{\partial r}+{\displaystyle \sum _{m=1}^M{\alpha }_m{p}_m}\text{on}{S}_R.$$E7

Denoting a normal-mode function at the radiation boundary by $$f_m(z)$$, which is associated with the m-th propagating mode eigenvalue, the pressure modes in Eq. (6) may be written as

$$p_m(z)=a_m{f}_m(z),m=1,2,\cdots ,M,$$E8

where $$a_m$$ denotes a modal participation factor. Then the radiation boundary condition in Eq. (7) may be rewritten as

$$p(z)+{\displaystyle \sum _{m=1}^M{a}_m{\alpha }_m{f}_m(z)}=0\text{on}{S}_R,$$E9

where the constants $$a_m$$ are determined by using the $$(1/\rho (z))$$ -orthogonality of the normal modes. It has tacitly been assumed here that the waveguide has constant water depth and range-independent but depth-dependent sound speed in the vicinity of the truncation boundary $$S_{\text{R}}$$ and beyond, so that the depth eigenproblem corresponding to the problem in Eq. (1) could be solved at least numerically [23].

Note that while the radiation condition in Eq. (4) on an individual mode is local, the radiation condition in Eq. (9) is global, meaning that nodes of an element on the truncation boundary are linked to other elements there in view of the coefficients $$a_m$$.

5. Constraints

In view of Eq. (8), Eq. (6) may be written as

$$[C]{\left\{p(z),a_1,a_2,\dots ,a_M\right\}}^T=0\text{on}{S}_R,$$E10

where

$$[C]=\left[1,-f_1(z),-f_2(z),\dots ,-f_M(z)\right],$$E11

and the companion vector in Eq. (10) is unknown. Equation (10) will be treated as a constraint in the following FE model.

6. Variational formulation

For the purpose of finite element modeling, it would be convenient to construct a variational formulation [34-35]. In the present study, in order to avoid possible numerical difficulties with handling a point source, a small fluid domain $${\mathrm{\Omega }}_S$$ surrounding the source has been excluded so that the computational domain is $${\overline{\mathrm{\Omega }}}_I={\mathrm{\Omega }}_I-{\mathrm{\Omega }}_S$$. Consider the following axisymmetric functional $$I(p)$$ defined in the cylindrical coordinate system $$(r,z)$$ (see Fig. 1):

$$I(p)=\frac{1}{2}{\displaystyle \underset{{\overline{\mathrm{\Omega }}}_I}{\int }\frac{1}{\rho }}\left\{{\left(\frac{\partial p}{\partial r}\right)}^2+{\left(\frac{\partial p}{\partial z}\right)}^2-k^2{p}^2\right\}r\text{}drdz+\frac{1}{2}{\displaystyle \sum _{m=1}^M{\displaystyle {\int }_{S_{\text{R}}}\frac{1}{\rho }{\alpha }_m{p}_m^2}}rdz-{\displaystyle \underset{S_{\text{D}}+S_{\text{N}}}{\int }\frac{1}{\rho }}\frac{\partial p}{\partial n}pdS,$$E12

where $$S_{\text{D}}$$ denotes the surface on which a Dirichlet boundary condition is prescribed and $$S_{\text{N}}$$ the surface with prescribed Neumann boundary condition, and the other domains of integration are identified in Fig. 1.

It can readily be shown that the variational condition

$$\delta I=0$$E13

leads to the governing differential equation in Eq. (1) and the boundary conditions in Eqs. (2)–(4). Thus, Eqs. (12) and (13) can be used to develop an FE model using the Rayleigh–Ritz approximation. However, the resulting solution should also obey the constraints in Eq. (10), which will ensure the imposition of the radiation boundary condition as discussed above. This will be achieved by modifying the discrete approximation to the functional in Eq. (12).

7. Finite element model

The finite fluid domain $${\overline{\mathrm{\Omega }}}_I$$ (which excludes the source) of the axisymmetric waveguide in Fig. 1 may be discretized using eight-noded axisymmetric quadrilateral elements with $$C_o$$ -continuity and the well-known isoparametric formulation [15]. The computational domain is discretized into a mesh of finite elements. The finite element approximation for the field variable p may then be written as

$$p(r,z)\approx {\displaystyle \sum _{j=1}^{\tilde{n}}{\overline{p}}_{ej}}{N}_j(\xi ,\zeta )={\left[N\right]}^T\left\{{\overline{p}}_e\right\},$$E14

where $$\tilde{n}$$ denotes the number of element nodes (eight in the present study), $${\overline{p}}_{ej}$$ the nodal pressure variable/degrees of freedom (dofs) and $$N_j(\xi ,\zeta )$$ the polynomial shape function in the parametric coordinates $$(\xi ,\zeta )$$ in the $$(r,z)$$ plane [15]. The subscript e is used to indicate the quantity at the element. Substituting Eq. (14) into Eq. (12) yields the following discrete form:

$$I(p_e)\approx \frac{1}{2}{\left\{{\overline{p}}_e\right\}}^T\left(\left[K_e\right]-\left[M_e\right]\right)\left\{{\overline{p}}_e\right\}+\frac{1}{2}{\displaystyle \sum _{m=1}^M{\left\{{\overline{p}}_{em}\right\}}^T\left[{\overline{R}}_{em}\right]}\left\{{\overline{p}}_{em}\right\}-{\left\{{\overline{p}}_e\right\}}^T\left\{f_e\right\},$$E15

where $$\left\{{\overline{p}}_{em}\right\}$$ denotes nodal pressure on the radiation boundary due to the m-th mode and the various matrices above will be identified subsequently. The stationary condition of the potential $$I(p_e)$$ above should be sought subject to the constraint in Eq. (10). There are two ways of implementing this, one the classical Lagrangian multiplier approach and the other the penalty function approach; the latter, which is commonly used in the context of finite element analysis [15, 13] is adopted in the present work. To achieve this, a modified potential $$I^{\prime }$$ may be defined as

$$I^{\prime }=I+\frac{1}{2}{\left\{{{\overline{p}}^{\prime }}_e\right\}}^T{\left[C_e\right]}^T\left[{\beta }_P\right]\left[C_e\right]\left\{{{\overline{p}}^{\prime }}_e\right\},$$E16

where $$\left[C_e\right]$$ denotes the constraint matrix in Eq. (11) specific to an element. The penalty coefficient matrix $$\left[{\beta }_P\right]$$ above may be chosen to be diagonal for convenience, with $${\beta }_{Pm}$$ denoting the penalty parameter associated with the m-th mode. Equation (16) may be expanded as

$$I^{\prime }({p^{\prime }}_e)=\frac{1}{2}{\left\{{{\overline{p}}^{\prime }}_e\right\}}^T\left(\left[{K^{\prime }}_e\right]-\left[{M^{\prime }}_e\right]+\left[{R^{\prime }}_e\right]+{\left[C_e\right]}^T\left[{\beta }_P\right]\left[C_e\right]\right)\left\{{{\overline{p}}^{\prime }}_e\right\}-{\left\{{{\overline{p}}^{\prime }}_e\right\}}^T\left\{{f^{\prime }}_e\right\},$$BB1

where the enlarged element dof vector is defined as

$$\left\{{{\overline{p}}^{\prime }}_e\right\}=\left\{\begin{array}{c}\left\{{\overline{p}}_e\right\}\\ \left\{a\right\}\end{array}\right\}.$$E17

The enlarged stiffness, mass and damping matrices, and load vector in expanded Eq. (16), consistent with $$\left\{{{\overline{p}}^{\prime }}_e\right\}$$ in Eq. (17), are given by

$$\left[{K^{\prime }}_e\right]=\left[\begin{array}{cc}\left[K_e\right]& 0\\ 0& 0\end{array}\right]\left[{M^{\prime }}_e\right]=\left[\begin{array}{cc}\left[M_e\right]& 0\\ 0& 0\end{array}\right]\left[{R^{\prime }}_e\right]=\left[\begin{array}{cc}0& 0\\ 0& \left[R_e\right]\end{array}\right]\left\{{f^{\prime }}_e\right\}=\left\{\begin{array}{c}\left\{f_e\right\}\\ 0\end{array}\right\}.$$E18

The matrices $$\left[K_e\right]$$, $$\left[M_e\right]$$, and $$\left[R_e\right]$$ in Eq. (18) are traditionally called the element stiffness, mass and radiation damping matrices, and $$\left\{f_e\right\}$$ the load vector, respectively. They are given as follows:

$$\begin{array}{cc}\left[K_e\right]={\displaystyle \underset{{\mathrm{\Omega }}_e}{\int }\frac{1}{\rho }}{\left[\nabla N\right]}^T\left[\nabla N\right]d\mathrm{\Omega }& (\text{a})\\ \left[M_e\right]={\displaystyle \underset{{\mathrm{\Omega }}_e}{\int }\frac{k^2}{\rho }}{\left[N\right]}^T\left[N\right]d\mathrm{\Omega }& (\text{b})\\ \left[R_e\right]=\text{diag}\left[R_{e1},R_{e2},\dots ,R_{eM}\right]& (\text{c)}\\ R_{em}={\left\{f_{zm}\right\}}^T\left[{\overline{R}}_{em}\right]\left\{f_{zm}\right\}\left[{\overline{R}}_{em}\right]={\alpha }_m{\displaystyle \underset{S_{Re}}{\int }\frac{1}{\rho }{\left[N\right]}^T\left[N\right]dS}& (\text{d)}\\ \left\{f_{zm}\right\}={⟨f_{zm}(z_1),f_{zm}(z_2),\dots ,f_{zm}(z_n)⟩}^T& (\text{e)}\\ \left\{f_e\right\}={\displaystyle \underset{S_{Ne}}{\int }\frac{1}{\rho }{p}_v{\left[N\right]}^{T}dS}& (\text{f)}\end{array}$$E19

where $$\left[N\right]$$ denotes the shape-function matrix [see Eq. (14)], $$p_{\nu }=\partial p/\partial n$$, and $$f_{zm}(z_j)$$ denotes the j-th nodal value of the m-th mode on a finite element in contact with the radiation boundary $$S_{Re}$$. The steps required to derive Eq. (19c) are outlined in Appendix A. It is of interest to note that the radiation-damping matrix $$\left[R_e\right]$$ in Eq. (19c) implies uncoupled modal participation. However, the constraint term involving the matrix $$\left[C_e\right]$$ in expanded Eq. (16) brings about modal coupling. The various integrals above are defined over relevant finite element domains. The stationary condition of the potential $$I^{\prime }$$ in expanded Eq. (16) is obtained by setting

$$\frac{\partial I^{\prime }}{\partial \left\{{{\overline{p}}^{\prime }}_e\right\}}=0$$E20

Equation (20) leads to general element equations of the form

$$\left(\left[{K^{\prime }}_e\right]-\left[{M^{\prime }}_e\right]+\left[{R^{\prime }}_e\right]+{\left[C_e\right]}^T\left[{\beta }_P\right]\left[C_e\right]\right)\left\{{{\overline{p}}^{\prime }}_e\right\}=\left\{{f^{\prime }}_e\right\}.$$E21

It may be noted that if the penalty matrix $$\left[{\beta }_P\right]=0$$ in Eq. (21), the constraints are ignored; as the penalty parameter values increase, the error in satisfying the constraint equations decreases, and for very high values of penalty, the numerical solution may break down. Hence, a judicious choice of the penalty parameters is essential. The radiation-damping matrix $$\left[{R^{\prime }}_e\right]$$ and the constraint matrix $$\left[C_e\right]$$ in Eq. (21) correspond to elements on the radiation boundary. Hence, for this case, the FE equation may be deduced from Eqs. (18) and (21) as

$$\left(\left[K_e\right]-\left[M_e\right]\right)\left\{{\overline{p}}_e\right\}=\left\{f_e\right\}.$$E22

The radiation-damping matrix $$\left[{R^{\prime }}_e\right]$$ in Eq. (19c), which is complex in view of Eq. (5), is defined only for elements that share one or more of their boundaries with the artificial boundary $$S_R$$. Carrying out thus the finite element assemblage operation, yields the following global finite element equations:

$$\left(\left[K^{\prime }\right]-\left[M^{\prime }\right]+\left[R^{\prime }\right]+\left[C^{\prime }\right]\right)\left\{{\overline{p}}^{\prime }\right\}=\left\{f^{\prime }\right\},$$E23

where the global solution vector $$\left\{{\overline{p}}^{\prime }\right\}$$ consists of all the pressure dof in the computational domain as well as the unknown vector $$\left\{a\right\}$$ in Eq. (8).

The global finite element matrices in Eq. (23) may formally be written as

$$\begin{array}{l}\left[K^{\prime }\right]={\displaystyle \sum _e\left[{K^{\prime }}_e\right]}\left[M^{\prime }\right]={\displaystyle \sum _e\left[{M^{\prime }}_e\right]}\left[R^{\prime }\right]={\displaystyle \sum _e\left[{R^{\prime }}_e\right]}\\ \left[C^{\prime }\right]={\displaystyle \sum _e{\left[C_e\right]}^T\left[{\beta }_P\right]\left[C_e\right]}\left\{f^{\prime }\right\}={\displaystyle \sum _e\left\{{f^{\prime }}_e\right\}}\end{array}$$E24

where $${\mathrm{\Sigma }}_e$$ denotes the standard finite element assemblage operation [15].

8. Modeling of a point source

When the inhomogeneous Helmholtz equation in Eq. (1) is employed in the FE model, the source term involving the delta function, as the other terms of the differential equation, is satisfied only approximately over the finite elements in contact with the point source. Of course, the error is expected to decrease with mesh refinement. The present FE formulation uses the complex pressure p as the field variable. Hence, a kinematic/Dirichlet boundary condition in terms of p would be satisfied exactly at the finite element nodes. In light of this, it would be interesting to see whether the effect of the source could be modeled as a kinematic boundary condition. To facilitate this, the computational domain employed above (see Eq. (12)) excludes the source. This is achieved by matching each finite element node with the source, and excluding all the finite elements that are in contact with the source node. Then the free field pressure due to the source on the periphery of the excluded domain is imposed as a kinematic boundary condition in the finite element model. The discontinuity of the fields on the periphery of the region enclosing the source is our equivalent source. It may be argued that the pressure distribution on the excluded domain boundary is not the actual one, which would be known only after solving the FE equations. However, the following argument justifies the approach. It is known that for small volume sources, the pressure in the far field is not affected by the individual shape of a source, as long as the source strengths are equal. Thus, this justifies the use of a computational domain that excludes a small FE domain around a point source. In the present study, the size of the excluded domain has been kept at about a tenth of the wavelength. Comparison of the FE results with an analytical solution indicates that such a choice is satisfactory.

9. Solution of FE equations

The global FE equation in Eq. (23) may be written for brevity as

$$\left[A\right]\left\{{\overline{p}}^{\prime }\right\}=\left\{f^{\prime }\right\}$$E25

It may be noted that for an acoustic medium with real sound speed, the coefficient matrix $$\left[A\right]$$ above is complex even though $$\left[K^{\prime }\right]$$, $$\left[M^{\prime }\right]$$, and $$\left[C^{\prime }\right]$$ are real. This is because $$\left\{f^{\prime }\right\}$$ is complex. Also, note that $$\left[R^{\prime }\right]$$ is complex due to the presence of $${\alpha }_m$$ in Eq. (19d). For a lossy medium modeled with complex sound speed, $$\left[M^{\prime }\right]$$ is also complex. Although $$\left[A\right]$$ is non-self-adjoint, it is a complex symmetric matrix and hence the Gauss solver employed here to obtain the solution to Eq. (25) exploits the attendant computational advantage. Since such solvers for FE equations are coded as block solvers with compact storage scheme, large finite element models can be handled even with modest computer storage. Of course, such a solution strategy involves overhead in the form of read/write operations on secondary storage devices. This approach may be contrasted against those of Bayliss et al. [36] and Athanassoulis et al. [37] who have used iterative methods based on the conjugate-gradient technique. Solvers based on the conjugate-gradient method have been found much more efficient than Gauss solvers when the size of the matrix equation is very large, say, several tens of thousands of equations, and hence they hold promise for high frequency FE models.

Since the present FE model adopts a penalty function approach to impose the radiation boundary condition with multiple radiating modes, the choice of suitable penalty parameter $${\alpha }_{Pm}$$ is important. This can be resolved through numerical experiments. The penalty parameter was obtained by prescribing a scale factor on the average value of the diagonals of the coefficient matrix $$\left[A\right]$$ in Eq. (25); i.e.,

$${\beta }_{Pm}=\frac{{\beta }_s}{n^{\prime }}{\displaystyle \sum _{i=1}^{n^{\prime }}\left|A_{ii}\right|},$$E26

where $$n^{\prime }$$ denotes the total number of FE equations/dof and $${\beta }_s$$ a user-specified penalty scale factor. Computations indicate that the results are stable over a wide range of $${\beta }_s$$ values. The results reported here have been obtained using $${\beta }_s=100.$$

10. Normal modes in an ocean waveguide with depth dependence

The sound speed in an ocean-acoustic waveguide is in general both depth- and range-dependent. Depth dependence is considered very important because it is responsible for many interesting phenomena in waveguide propagation. The two well-known methods that have been developed to study acoustic waves in depth-dependent waveguides are the fast-field technique and the normal-mode expansion [9, 1], the latter being the method that we have used. The normal-mode approach consists of first solving the depth eigenproblem for a given sound-speed profile to obtain the radial wavenumbers and the associated depth modes, which respectively are the eigenvalues and eigenfunctions. The depth eigenproblem could be solved exactly only for a few special profiles. In the finite element model for depth- and range-dependent waveguides [31-32], the eigensolution of the depth problem is required for imposition of the radiation condition at the truncation boundary. For such applications, it would be convenient to have the depth modes in a compact analytical form. We have explored this aspect with specific reference to shallow-water waveguides.

The depth eigenproblem can be cast in a variational form by suitably defining a functional. Then, the classical Rayleigh–Ritz method may be employed to find a variational approximation to the eigensolution of the depth problem in ocean-acoustic waveguides. The depth modes obtained would have a more compact analytical form than those derived using finite difference or finite element methods. The present work provides an RR model for the depth eigenproblem and demonstrates its utility for shallow-water waveguides.

11. Mathematical model

For the cylindrically symmetric waveguide having depth-dependent density $$\rho$$ and sound speed $$c$$, the inhomogeneous pseudo Helmholtz equation governing the linear harmonic acoustic pressure field $$p(r,z)$$ in the waveguide is given in cylindrical coordinates (r,z) as [9, 1]

$$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial p}{\partial r}\right)+\rho (z)\frac{\partial }{\partial z}\left(\frac{1}{\rho (z)}\frac{\partial p}{\partial z}\right)+\frac{{\omega }^2}{c^2(z)}p=-\frac{\delta (r)\delta (z-z_s)}{2\pi r},$$E27

where r denotes the range coordinate and z the depth coordinate as shown in Fig. 2, and the r.h.s. denotes a point source of unit amplitude located at r = 0 and $$z=z_s$$, with $$\delta$$ denoting the Dirac delta function. Eq. (27) can also be applied to problems with attenuation by introducing a complex sound speed.

Figure 2.

A two-layer cylindrically symmetric waveguide

A variable separable solution for the homogeneous form of Eq. (27) may be written as

$$p(r,z)=\overline{R}(r)Z(z)$$E28

Then, upon using Eq. (28) in the homogeneous form of Eq. (27), the following ordinary differential equations are obtained:

$$\frac{1}{r}\frac{d}{dr}\left(r\frac{d\overline{R}(r)}{dr}\right)+k_r^2\overline{R}(r)=0$$E29

$$\rho (z)\frac{d}{dz}\left(\frac{1}{\rho (z)}\frac{dZ(z)}{dz}\right)+\left(\frac{{\omega }^2}{c^2(z)}-k_r^2\right)Z(z)=0$$E30

where $$\omega$$ denotes the circular frequency and $$k_r^2$$, the separation constant, which turns out to be the square of the radial/horizontal wavenumber. Eq. (29) evidently pertains to the radial/horizontal modes $$\overline{R}\left(r\right)$$, and Eq. (30) pertains to the depth modes $$Z\left(z\right)$$. Choosing a pressure release boundary at the top $$\left(z=0\right)$$ and a mixed/Robin boundary condition at the seabed $$\left(z=D_1\right)$$, the boundary conditions of our problem are written as [1, 30, 33]

$$Z(0)=0$$E31

$$Z(D_1)+\frac{g\left(k_r^2\right)}{\rho }\frac{dZ(D_1)}{dz}=0$$E32

where

$$g(k_r^2)={\rho }_b/\sqrt{\left(k_r^2-\frac{{\omega }^2}{c^2}\right)}$$BB2

with $${\rho }_b$$ denoting the density of the acoustic fluid in the isovelocity half-space underlying the water column. Eq. (32) facilitates replacing the half-space in the Pekeris waveguide [38] by means of an impedance-type boundary condition. It may be noted that Eq. (30) together with the homogeneous boundary conditions in Eqs. (31) and (32) do not constitute a proper Sturm-Liouville problem because Eq. (32) depends on the unknown eigenvalue $$k_r^2$$. Porter and Reiss [30] employed a finite difference model to solve Eq. (30) together with the boundary conditions in Eqs. (31) and (32). As an alternative to the above formulation, the waves in the fluid half-space below the water column are also considered here [1]. The governing equation for the waves in this fluid half-space is given as,

$${\rho }_{\text{b}}(z)\frac{d}{dz}\left(\frac{1}{{\rho }_{\text{b}}(z)}\frac{d{Z}_{\text{b}}(z)}{dz}\right)+\left(\frac{{\omega }^2}{c_{\text{b}}^2(z)}-k_r^2\right)Z_{\text{b}}(z)=0,\text{}{D}_1\le z\le \infty$$E33

where $$Z_{\text{b}}$$ denotes the depth function in the fluid half-space having depth-dependent density $${\rho }_{\text{b}}$$ and sound speed $$c_{\text{b}}$$. The interface conditions at the seabed are given by the kinematic and dynamic conditions,

$$\begin{array}{l}Z(D_1)=Z_{\text{b}}(D_1)(\text{a})\\ \frac{1}{\rho }\frac{dZ(D_1)}{dz}=\frac{1}{{\rho }_b}\frac{d{Z}_b(D_1)}{dz}\text{(b})\end{array}$$E34

In addition, the depth mode $$Z_{\text{b}}$$ should remain bounded as $$z\to \infty$$. Our primary objective is to consider a variational formulation for Eqs. (30) and (33), together with appropriate boundary conditions, and obtain a RR approximation to the depth-dependent problem.

12. Variational formulation and Rayleigh–Ritz approximation

A variational formulation that leads to the boundary value problem in the last section is sought now. The operator being symmetric, there exists a functional, the variation of which leads to Eq. (30) and appropriate boundary conditions, and similarly for the half-space. Consider the functional $$\Pi (Z)$$ and $${\Pi }_b(Z_b)$$ defined respectively in the water column and the half-space as

$$\prod (Z)=\frac{1}{2}{\displaystyle \underset{0}{\overset{D{}_1}{\int }}\left[\frac{1}{\rho (z)}{\left(\frac{dZ}{dz}\right)}^2-\frac{{\omega }^2}{\rho (z)c^2(z)}{Z}^2+\frac{1}{\rho (z)}{k}_r^2(z)Z^2\right]}dz{-\frac{1}{\rho (z)}{Z}_{\upsilon }Z|}_0^{D_1}$$E35

$$\begin{array}{l}{\prod }_b(Z_b)=\frac{1}{2}{\displaystyle \underset{D_1}{\overset{D_2}{\int }}\left[\frac{1}{{\rho }_b(z)}{\left(\frac{d{Z}_b}{dz}\right)}^2-\frac{{\omega }^2}{{\rho }_b(z)c^2(z)}{Z}_b^2+\frac{1}{{\rho }_b(z)}{k}_r^2{Z}_b^2\right]dz}\\ -\frac{1}{{\rho }_b(z)}{Z}_{b\upsilon }(D_1)Z_b(D_1)+\frac{1}{2}\beta \frac{Z_b^2(D_2)}{{\rho }_b},D_2\to \infty \end{array}$$E36

where suffix $$\nu$$ denotes z-derivative. At the interface $$z=D_1$$ between the water column and the half-space, the conditions noted in Eq. (30) must be imposed. In view of Eq. (34b), this can be achieved by setting in Eq. (35)

$$\frac{1}{\rho }{Z}_{\nu }(D_1)=\frac{1}{{\rho }_b}\frac{d{Z}_b(D_1)}{dz},$$E37

and in Eq. (36),

$$\frac{1}{{\rho }_b}{Z}_{b\nu }(D_1)=\frac{1}{\rho }\frac{dZ(D_1)}{dz}$$E38

In addition, Eq. (34a) should be imposed. Then, it can be easily shown that the variational condition $$\delta \Pi =0$$ leads to Eq. (30) and the boundary conditions in Eq. (31) as well as the interface condition in Eq. (37), where $$\delta$$ denotes the first variation. Similarly, the variational condition $$\delta {\Pi }_b=0$$ leads to Eq. (33), and the interface conditions in Eq. (34a) and Eq. (38). In addition, at $$z=D_2$$, we obtain the condition

$$\frac{d{Z}_b(D_2)}{dz}+\beta Z_b(D_2)=0,$$E39

where (see Eq. (32)),

$$\beta =\sqrt{\left(k_r^2-\frac{{\omega }^2}{c_b^2}\right)}$$BB3

Note that there are three cases that can be analyzed using Eq. (36):

Case 1: $$D_2$$ is finite, $$\beta \to 0$$

This corresponds to the case when a depth-dependent seabed of finite thickness is terminated by a rigid boundary.

Case 2: $$\beta$$ is finite, $$D_2\to \infty$$

This corresponds to a depth-dependent seabed of infinite thickness.

Case 3: $$D_2$$ and $$\beta$$ are finite

This is a three-layer problem, where the top layer is the water column, the second layer is a layer of seabed with depth varying density and speed, and the bottom layer represents the seabed of infinite extent with uniform sound speed and density.

We now seek an assumed mode solution with n terms to the above variational problem in the form

$$\begin{array}{ccc}Z(z)\approx {\displaystyle \sum _{j=1}^{n_1}{\overline{\varphi }}_{{}_j}{\psi }_j(z)}& 0\le z\le D_1& (\text{a})\\ Z_b(z)\approx {\displaystyle \sum _{j=1}^{n_2}{\overline{\varphi }}_{bj}{\psi }_{bj}(z)}& D_1\le z\le \infty & (\text{b})\end{array}$$E40

where $$n=n_1+n_2$$, and $${\psi }_j$$ and $${\psi }_{bj}$$ denote the known mode function (coordinate function) satisfying the kinematic boundary condition in the water column and $${\overline{\varphi }}_j$$ an unknown constant, and their counterparts with suffix b correspond to those of the half-space. The two sets of mode functions above are such that they satisfy the relevant boundary conditions as well as the interface conditions in Eq. (34) and the conditions in Eq. (39). Such functions may readily be constructed by solving a two-layer depth problem, which is nothing but the Pekeris waveguide [38], with an appropriate choice of constant velocity and density in the water column and the seabed half-space. This approach has been adopted here. Then, it follows that the continuity of pressure field at the interface $$z=D_1$$ implies that the assumed mode expansion in Eq. (40) reduces as,

$$Z(z)\approx {\displaystyle \sum _{j=1}^n{\overline{\varphi }}_{{}_j}{\psi }_j(z)}0\le z\le D_2$$E41

where we have combined the depth modes of a two-layer isovelocity waveguide as one combined set with redefined coefficients $$\overline{\varphi }$$. Then, using Eq. (40) in the functionals in Eqs. (35) and (36), and combining them, an algebraic approximation for the functionals is obtained as

$$\overline{\prod }=\prod +{\prod }_b=\frac{1}{2}{\left\{\overline{f}\right\}}^T\left[K^I\right]\left\{\overline{f}\right\}-\frac{1}{2}{\omega }^2{\left\{\overline{f}\right\}}^T\left[K^{II}\right]\left\{\overline{f}\right\}+\frac{1}{2}{k}_r^2{\left\{\overline{f}\right\}}^T\left[K^{III}\right]\left\{\overline{f}\right\}$$E42

where

$$\begin{array}{l}{K}_{ij}^I={\displaystyle \underset{0}{\overset{D_1}{\int }}\frac{1}{\rho (z)}\frac{d{\psi }_i}{dz}}\frac{d{\psi }_j}{dz}dz+{\displaystyle \underset{D_1}{\overset{D_2}{\int }}\frac{1}{{\rho }_b(z)}\frac{d{\psi }_{bi}}{dz}}\frac{d{\psi }_{bj}}{dz}dz\\ K_{ij}^{II}={\displaystyle \underset{0}{\overset{D_1}{\int }}\frac{1}{\rho (z)c^2(z)}{\psi }_i{\psi }_j}dz+{\displaystyle \underset{D_1}{\overset{D_2}{\int }}\frac{1}{{\rho }_b(z)c_b^2(z)}{\psi }_{bi}{\psi }_{bj}}dz\\ K_{ij}^{III}={\displaystyle \underset{0}{\overset{D_1}{\int }}\frac{1}{\rho (z)}{\psi }_i{\psi }_j}dz+{\displaystyle \underset{D_1}{\overset{D_2}{\int }}\frac{1}{{\rho }_b(z)}{\psi }_{bi}{\psi }_{bj}}dz\end{array}$$E43

It has been assumed in the above that the contribution due to the term in Eq. (36) is negligible as $$D_2\to \infty$$. Further, it may be noted that since the boundary and interface conditions are satisfied by the trial functions chosen above, when the functionals in Eqs. (35) and (36) are combined to obtain Eq. (42), the boundary and interface terms add up to become trivial and hence do not contribute to the discrete approximation in Eq. (42).

The variational condition is now replaced by the condition

$$\frac{\partial \overline{\prod }}{\partial {\varphi }_j}=0,j=1,2,\mathrm{...},n$$E44

Eq. (44) yields a symmetric algebraic eigenproblem given by

$$\left({\omega }^2\left[{\text{K}}^{\text{II}}\right]-\left[{\text{K}}^{\text{I}}\right]\right)\left\{\overline{\varphi }\right\}=k_r^2\left[{\text{K}}^{\text{III}}\right]\left\{\overline{\varphi }\right\}$$E45

The eigensolution of Eq. (45) may be denoted as

$$\left(k_{rj}^2,\left\{{\varphi }^{\left(j\right)}\right\}\right),j=1,2,\mathrm{..}n$$E46

It may be noted here that the eigenproblem in Eq. (45) remains linear unlike the Porter and Reiss model that is based on Eqs. (30)–(32). Having obtained the eigenvalues $$k_{rj}^2$$ and the eigenvectors $$\left\{{\overline{\varphi }}^{(j)}\right\}$$, the eigenfunctions /depth modes may be written, using Eq. (40), as

$${{\rm Z}}_j\left(z\right)={\left\{{\overline{\varphi }}^{(j)}\right\}}^T\left\{\psi \left(z\right)\right\}$$E47

where

$${\left\{\psi \left(z\right)\right\}}^T=⟨{\psi }_1(z),{\psi }_2(z),\cdots ,{\psi }_n(z)⟩$$E48

Eq. (47) provides a compact semi-analytical form for the depth modes that are convenient to employ in FE models such as those in Fix and Marin [31] and Vendhan et al. [32] for approximating the radiation condition at the truncation boundary. The depth modes obtained can of course be used to set up the normal-mode solution to the forced Helmholtz equation in Eq. (27). Note that for a Pekeris waveguide, the normal-mode solution based on the discrete spectrum has to be augmented with the continuous spectrum contribution [1]. Since the eigenvectors in Eq. (45) are $$\left[{\text{K}}^{\text{III}}\right]$$ -orthogonal, it can easily be shown that the eigenfunctions in Eq. (47) satisfy the following orthogonality condition:

$${\displaystyle \underset{0}{\overset{D_1}{\int }}\frac{1}{\rho (z)}{Z}_i{Z}_j}dz+{\displaystyle \underset{D_1}{\overset{D_2}{\int }}\frac{1}{{\rho }_b(z)}{Z}_i{Z}_j}dz=0,i\ne j$$E49

The orthonormal depth functions are obtained as

$${\overline{Z}}_j(z)=Z_j/\sqrt{{\left\{{\overline{\varphi }}^{(j)}\right\}}^T\left[K^{III}\right]\left\{{\overline{\varphi }}^{(j)}\right\}}$$E50

In terms of finite element terminology, the RR model for each layer may be looked upon as a super-element with $$C_1$$ continuity at the inter-element boundary and the operation leading to Eq. (42) is equivalent to element-assemblage operation.

13. Numerical analysis and discussion

In our Rayleigh–Ritz model, the first task is to compute the symmetric matrices $$\left[{\text{K}}^{\text{I}}\right]$$, $$\left[{\text{K}}^{\text{II}}\right]$$, and $$\left[{\text{K}}^{\text{III}}\right]$$ in Eq. (43). The next task is to find the eigensolution to Eq. (45). For problems with no attenuation, the real eigenvalues have been obtained employing the bisection method. For problems with attenuation, approximations to the complex roots have been obtained using a search procedure [39] and the eigenvalues refined by employing Newton–Raphson iteration. In all cases, the eigenvectors are obtained using inverse iteration.

To validate our algorithm, we applied the Rayleigh–Ritz model first to single-layer isovelocity waveguide examples without attenuation for which exact solutions are available. Different sound-speed profiles have been chosen to evaluate the accuracy of the RR model. Attenuation in the fluid half-space has also been considered. Different sets of RR approximations have been obtained by varying the number of assumed modes n in Eq. (41). The results for n = 2n_p, where n_p denotes the number of propagating modes turned out to be of good accuracy.

One should note the following remarks in connection with the performance of the RR model for the depth eigenproblem:

1. The mode shapes of an isovelocity waveguide have been chosen as trial functions, which satisfy appropriate interface conditions and the condition at the free surface. This renders the RR matrix highly diagonally dominant, which also greatly aids in numerical evaluation of the eigensolution.

2. For ocean waveguides, the depth variation of the sound speed is normally only a small percentage of the unperturbed value.

3. When the variation in sound speed is large, the above procedure may not give good results. One has to resort to high-order solutions. Even then, one can expect accurate eigenvalues, but not eigenvectors. This is because the convergence rate for the eigenvectors is slower than that for the eigenvalues.

14. Numerical examples

We considered several examples to illustrate the versatility of our FEM approach in remote sensing problems. In all our examples, we employed a Dirichlet boundary condition on the air–sea interface, a Neumann boundary condition on the ocean bottom boundary, and a unit point source at a depth of 36 m below the air–water interface. Both depth-dependent and uniform sound-speed water columns are considered.

14.1. Isovelocity case

The finite element method for the solution of inhomogeneous ocean-acoustic waveguide problems is validated first with analytical results for isovelocity waveguides. A cylindrically symmetric plane parallel waveguide of depth 100 m with a point source is shown in Fig. 3. The finite element model consists of a uniform grid of isoparametric quadrilateral elements, with the element length being about a tenth of the source signal wavelength. As discussed previously, a domain of two elements has been excluded to remove the source from the truncated domain (Fig. 1). The FE mesh consists of 1000 elements in range and 60 elements in depth. The computed acoustic pressure along the range at the depth of the source is compared in Fig. 3 with the normal-mode solution with 50 modes, of which only the leading few modes are propagating. In all cases, the FEM results compared well with analytical results. The mesh is chosen appropriately so that the modal error is less than 5%.

Figure 3.

Idealized ocean waveguide

14.2. Rectangular hump

Sea mounts are often encountered in under-water ocean problems. In order to understand their impact on wave propagation characteristics in shallow-water environment, we considered a rigid rectangular hump of width 40 m and height 20 m on ocean bottom as shown in Fig. 4(a). The contour map of transmission loss (TL) of a 60-Hz point source located on the z-axis at a depth of 36 m from the water surface is shown in Fig. 4. Panel (a) shows the TL in the presence of a rectangular hump on the seabed. Panel (b) shows the TL of the water column without the rectangular hump. Notice that the rectangular hump has a distinct signature in TL pattern especially on the right of the hump.

Figure 4.

Transmission loss of a shallow-water column with a (a) rectangular hump on the seabed and (b) flat-bottom surface

It is instructive to take a look at the acoustic power distribution in the modes. Figure 5 shows the modal power spectrum of the shallow-water column with the rectangular hump in panel (a) and without the rectangular hump in panel (b). Notice that there is a substantial redistribution of power among the modes due to the presence of the rectangular hump.

Figure 5.

Spectra of modal efficiencies for (a) seabed with a rectangular bump and (b) a flat seabed

14.3. Down-sloping bottom

Shallow-water conditions are encountered in the near-coast context where the ocean bottom has a sloping geometry. There are two situations to consider, up-slope and down-slope, depending on the location of the source with respect to the slope. First we consider the down-sloping case where the ocean-bottom slopes down from 100 m to 230 m over a distance of 600 m. The details of the geometry are shown in Fig. 6.

Figure 6.

Transmission loss of a shallow-water column with (a) down-sloping bottom and (b) flat bottom of depth 100 m

Figure 7.

Comparison of TL at 36 m depth of an ocean with down-sloping bottom and a flat bottom

Panel (a) shows the TL with the down-sloping bottom. Panel (b) shows the TL for a water column with the flat bottom at depth 100 m. Both results are for the source frequency of 150 Hz. Notice the distinct spatial power distribution manifested by the sloping bottom. To facilitate a better comparison, we have shown in Fig. 7 the TL at 36 m depth corresponding to the flat and sloping bottoms. Notice that the TLs for the two cases are similar in the region between the source and the middle of the slope. Beyond that, the TL corresponding to the sloping bottom is significantly larger than that of the flat bottom.

Figure 8.

Modal power spectrum of shallow-water column with (a) down-sloping bottom and (b) flat bottom of depth 100 m

In order to better understand the propagation phenomenology, the modal power spectrum for the shallow-water ocean with (a) down-sloping bottom and (b) flat bottom are shown in Fig. 8. Notice that there is a significant redistribution of energy in the case of sloping bottom although the total power flows in both cases are approximately the same.

14.4. Up-sloping bottom

Next we consider the problem of sloping bottom where the ocean bottom slopes up (with respect to location of the source) from 230 m to 100 m over a distance of 600 m. The details of the geometry are shown in Fig. 9. The acoustic source is located at 36 m below the water surface on the left.

Figure 9.

Transmission loss of shallow-water column with up-slope bottom at 105 Hz and 150 Hz

Panel (a) shows TL for the case of 105Hz and panel (b) shows the case of 150 Hz. We notice that at 105 Hz there is a substantial reduction in power flow. However, at 150 Hz the power flow is as good as that of a flat-bottom waveguide. The mode spectral distribution in Fig. 10 shows the details of how the power flows in the two cases. We notice that for the up-slope case, power flow can be good at certain frequencies and not good at others, depending on the impedance matching conditions. In contrast, for the case of down slope the power flow is good for all the frequencies that we studied.

Figure 10.

Modal spectrum of shallow-water column with up-sloping bottom (a) 105 Hz (b) 150 Hz

14.5. Object in the water column

Characterizing the signatures of objects in the ocean is an important remote sensing problem. We consider a cylindrical rigid object of radius 20 m in the middle of a water column as shown in Fig. 11. Panel (a) shows the TL for the source frequency at 135 Hz. Panel (b) shows the TL for the source frequency at 150 Hz. We notice that the power flow can be substantially influenced by the object, depending on the frequency of operation. This is because of the interference phenomena involving the object and the boundaries of the waveguide.

Figure 11.

TL of a shallow-water column with an object (a) 135 Hz (b) 150 Hz

14.6. Shallow-water column with rippled top surface

Ripples on the water surface can be generated by gravity and wind conditions. Such surface undulations can considerably influence the wave propagation in the shallow-water waveguide. To illustrate this phenomenon, we have taken a periodic structure on the air–water interface as shown in Fig. 12. The top surface has a sinusoidal undulation of amplitude 5 m and period 50 m.

Figure 12.

TL in a shallow-water column with (a) wind-generated rippled air–water interface, and (b) flat water surface

Our FEM results show that the surface ripples causes substantial transmission loss compared to that of the flat water surface for the case when the source frequency is 120 Hz. However, this pattern is quite sensitive to source frequency. For some frequencies, the TL may be large and for others, TL can be low. Dimensions of the waveguide and ripple geometry in terms of the source signal wavelength are key factors influencing the physics.

14.7. Shallow-water column with depth-dependent sound speed

In all the examples considered thus far, we have assumed that the water column has uniform sound speed. This is rarely true in practice even for the shallow-water ocean. The normal modes for the depth-dependent waveguide are required to impose the radiation boundary condition in our finite element procedure. The Rayleigh–Ritz approximation is used for obtaining normal modes for this problem. The sound-speed profile taken for this study is shown in Fig. 13.

Figure 13.

Sound-speed profile (S1) used for our study

Figure 14.

TL for the shallow-water ocean with (a) depth-dependent sound speed, and (b) uniform sound speed

The TL for our geometry with the sound-speed profile given in Fig. 13 is shown in Fig. 14. The result for source frequency of 150 Hz is shown in Panel (a) and that corresponding to isovelocity is shown in Panel (b). Although the sound-speed variation is very small, we notice the impact of depth dependence of sound speed on TL is substantial. However, at lower frequencies, this kind of sound-speed variation does not influence the TL much.

14.8. Shallow-water column with depth-dependent sound speed and a rectangular bump on seabed

Next, we consider the case of shallow-water ocean with depth-dependent sound speed and a rigid rectangular hump on the seabed.

Figure 15.

TL for the shallow-water ocean with (a) depth-dependent sound speed and a rectangular bump on the seabed, and (b) depth-dependent sound speed and flat bottom

We observe that, for the case when the source frequency is 60Hz, the presence of the small rectangular hump on the seabed has a significant effect on the transmission loss of a depth-dependent ocean.

14.9. Shallow-water column with depth-dependent sound speed and rippled top surface

Finally, we consider the case of shallow-water ocean with depth-dependent sound speed (Fig. 13) and a rippled air–water surface. The results are shown in Fig. 16.

Figure 16.

TL for the shallow-water ocean with (a) depth-dependent sound speed and a rippled air–water interface, and (b) depth-dependent sound speed

Note that for source frequency of 30Hz, the presence of ripples has reduced the transmission loss in most regions. This is in contrast to the last case (Fig. 15) where there is a rectangular bump on the seabed. However, these characteristics are due to interference phenomenon and hence have strong frequency dependence. The important point is that small features such as ripples can have a significant impact on the underwater propagation characteristics.

15. Conclusion

A finite element approach has been presented for remote sensing in shallow-water ocean environment. The three principal elements of remote sensing are: (a) signal propagation and reception, (b) data analysis, and (c) inversion or retrieval. This chapter exclusively deals with part (a) of the trilogy of remote sensing. Although several approaches have been developed for wave propagation studies in underwater ocean, they all have limitations when encountered with complex geometries and environments as in shallow-water ocean. An FE approach is both accurate and feasible for such applications. In order to minimize the problem size, a Bayliss-type damper was imposed to truncate the solution domain. Since several propagating modes can exist in the ocean waveguide, a penalty function approach was used to impose the radiation boundary condition in the variational finite element formulation of the problem. This penalty function approach was found to be robust over a wide range of penalty scale factors.

For the shallow-water ocean waveguide with depth-dependent sound-speed problem, the eigensolution was obtained using a Rayleigh–Ritz approximation. The trial functions are derived from an isovelocity problem that has exact solution. It is important to note that such trial functions automatically satisfy even the dynamic interface condition at the seabed, thus contributing to the accuracy of the numerical model. The proposed model is accurate and provides a compact semi-analytical form for the depth modes.

We thus have an accurate FE model for the remote sensing in range- and depth-dependent ocean-acoustic waveguides. Numerous examples were considered to illustrate the accuracy and versatility of this model. Admittedly, the computational effort in setting up the matrix in the proposed RR model using numerical quadrature is high compared to setting up the finite-difference-based matrix in the Porter and Reiss approach. However, noting the diagonal dominance of the matrix obtained in the RR model, it would be worthwhile exploring the possibility of approximating it by a narrow banded matrix in order to reduce the volume of computation in setting up the matrix and possibly in obtaining the eigensolution. We have also extended this work for the case of irregular elastic seabed. We continue to explore and further develop our finite element approach by applying it to several other ocean-acoustic problems encountered in the remote sensing of ocean environment.

16. Appendix A: Derivation of multimode radiation damping matrix

Consider the functional in Eq. (12). The contribution, $$I_R(p_e)$$, from the radiation boundary of a finite element is represented by the second integral in that equation; i.e.,

$$I_R(p_e)=\frac{1}{2}{\displaystyle \sum _{m=1}^M{\displaystyle \underset{S_{Re}}{\int }\frac{1}{\rho }}}{\alpha }_m{p}_m^{2}rdz$$E51

where M denotes the number of propagating modes, $${\alpha }_m$$ the damper coefficient associated with the m-th mode [see Eq. (5)], $$p_m(z)$$ the pressure associated with the m-th normal mode, and $$S_{Re}$$ the element surface on the radiation boundary (see Fig. 1).

The modal pressure on the radiation boundary is given by Eq. (8):

$$p_m(z)=a_m{f}_m(z)m=1,2,\dots ,M,$$E52

where $$f_m(z)$$ denotes the normal-mode function and $$a_m$$ the modal coefficient. Using the finite element representation, the modal pressure on the radiation boundary may be written as

$$p_m(z)={\displaystyle \sum _e{\left[N\right]}^T}\left\{{\overline{p}}_{em}\right\},$$E53

where $$\left[N\right]$$ denotes the shape functions and $$\left\{{\overline{p}}_{em}\right\}$$ the nodal pressure vector on an element edge on the radiation boundary due to the m-th mode. The summation symbol is used to indicate that Eq. (53) is a piecewise polynomial representation over the entire depth of the waveguide. Using Eqs. (52) and (53), Eq. (51) may be written in a discrete form for a finite element as [also see Eq. (15)]

$$I_R(p_e)=\frac{1}{2}{\displaystyle \sum _{m=1}^M{\left\{{\overline{p}}_{em}\right\}}^T\left\{{\overline{R}}_{em}\right\}\left\{p_{em}\right\}},$$E54

where

$$\left[{\overline{R}}_{em}\right]={\alpha }_m{\displaystyle \underset{S_{Re}}{\int }\frac{1}{\rho }}{\left[N\right]}^T\left[N\right]dS.$$E55

In view of Eq. (52), the vector of modal pressure at the nodes of an element in Eq. (53) may be written as

$$\left\{{\overline{p}}_{em}\right\}=a_m{⟨f_{zm}(z_1),f_{zm}(z_2),\dots ,f_{zm}(z_n)⟩}^T=a_m\left\{f_{zm}\right\},$$E56

where $$f_{zm}(z_j)$$ denotes the j-th nodal value of the m-th eigenmode on a finite element in contact with the radiation boundary. Now, using Eqs. (55) and (56), the functional in Eq. (54) may be written as

$$I_R=\frac{1}{2}{\displaystyle \sum _{m=1}^M{a}_m^2{R}_{em}},$$E57

where

$$R_{em}={\left\{f_{zm}\right\}}^T\left[{\overline{R}}_{em}\right]\left\{f_{zm}\right\}.$$E58

The foregoing steps form the basis for Eq. (19c).

17. Appendix B: Normal-mode functions for isovelocity waveguides

The Rayleigh–Ritz model presented for the depth eigenproblems employs the analytical depth modes for an isovelocity waveguide as the trial functions. The details of the various isovelocity waveguide examples encountered in the ocean context are presented here. It should be kept in mind that in our problem, the acoustic source and reception points are both in the water column. Therefore, the wave functions given here have been chosen particularly for this application.

For a single-layer waveguide of depth D with Dirichlet boundary condition on top and Neumann boundary condition on the bottom surface, the trial functions are given by

$${\psi }_j=a_j\sin (k_{zj}z)$$E59

where

$$k_{zj}=\frac{\left(j-0.5\right)\pi }{D},j=1,2,\mathrm{...}$$E60

where $$a_j$$ is chosen to normalize the mode functions.

For a two-layer waveguide shown in Fig. 17, the trial functions are given by

$$\begin{array}{c}{\psi }_j(z)=a_j\sin k_{zj}z0\le z\le D_1\\ =a_j\frac{\sin k_{zj}{D}_1}{\cos k_{bzj}(D_2-D_1)}\cos k_{bzj}(D_2-z)D_1\le z\le D_2\end{array}$$E61

where $$k_{zj}^2={\omega }^2/c_{}^2-k_{rj}^2,k_{bzj}^2={\omega }^2/c_b^2-k_{rj}^2$$ and $$k_{rj}$$ is the solution of the transcendental equation

$$\frac{k_{zj}{\rho }_b}{k_{bzj}\rho }=\tan k_{zj}{D}_1\tan k_{bzj}(D_2-D_1)$$E62

Figure 17.

A bounded two-layer waveguide

In the case of a Pekeris waveguide (for which $$D_2\to \infty$$ in Fig. 17), the trial functions are given by

$$\begin{array}{c}{\psi }_j(z)=a_j\sin k_{zj}z0\le z\le D_1\\ =a_j\sin (k_{bzj}{D}_1)e^{-i{k}_{bzj}(D_1-z)}{D}_1\le z\le \infty \end{array}$$E63

where $$k_{zj}^{}=\sqrt{{\omega }^2/c_{}^2-k_{rj}^2},k_{bzj}^{}=i\sqrt{k_{rj}^2-{\omega }^2/c_b^2}$$ and $$k_{rj}$$ is the solution of the transcendental equation

$$\tan k_{zj}{D}_1=-\frac{i{\rho }_b{k}_{zj}}{\rho k_{bzj}}$$E64

Note that discrete guided modes in the water column exist only for the case when $$k_{rj}^2>{\omega }^2/c_b^2$$.

There are two cases to consider:

Case 1: $$k<k_b$$

This applies to the situation when the sound-speed velocity in the seabed is smaller than that in the water column. In this case, there are no guided modes. The entire spectrum is continuous and does not have contribution to sound transmission in the water column at long distances.

Case 2: $$k>k_b$$

This applies to the situation when the sound speed in the seabed is larger than that in the water column. Here the spectrum consists of (a) discrete guided modes, (b) continuous radiation modes, and (c) surface modes. Among the three, it is the discrete guided modes that carry the sound signal over long distances in the water column.

Since our interest is in long-range sound transmission in the water column, we have restricted attention to discrete guided modes as shown above.

One should observe that our two-layer waveguide problem does not share the above mentioned behavior. Note that the two-layer waveguide is terminated at the bottom by a rigid boundary. Therefore, the underlying physical processes are different.

Case 1: $$k<k_b$$

Here the entire spectrum in the waveguide consists of discrete guided modes.

Case 2: $$k>k_b$$

In this case, the spectrum consists of discrete guided modes and surface modes. However, for long range propagation, the modes of significance are the discrete guided modes.

Acknowledgments

S. Mudaliar thanks the AFOSR for support. C.P. Vendhan thanks AOARD for support during this work period. The authors thank A.D. Chowdhury for his help in eigenvalue computations for depth-dependent water column.