**Remote Sensing of the Ocean Environment Using Finite Element Methods**

Saba Mudaliar, C.P. Vendhan and C. Prabavathi

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/62193

#### **Abstract**

[122] Lyzenga DR: Remote sensing of bottom reflectance and water attenuation parame‐ ters in shallow water using aircraft and Landsat data. Int J Remote Sens. 1981;2(1):

[123] Lyzenga DR: Shallow-water bathymetry using combined lidar and passive multi‐

[124] Minghelli-Roman A, Goreac A, Mathieu, S, Spigai M, Gouton P: Comparison of bathymetric estimation using different satellite images in coastal sea waters. Int J Re‐

[125] Teodoro AC, Gonçalves H, Pais-Barbosa J: Bathymetric estimation through principal components analysis using IKONOS-2 data, in Remote Sensing for Agriculture, Eco‐ systems, and Hydrology XII, Christopher M. U. Neale; Antonino Maltese, Editors,

[126] Eugenio F, Marcello J, Javier Martin J: High-resolution maps of bathymetry and benthic habitats in shallow-water environments using multispectral remote sensing

Proceedings of SPIE Vol. 7824 (SPIE, Bellingham, WA 2010), 782419.

imagery. IEEE Trans Geosci Remote Sens. 2015;53(7):3539-3549.

spectral scanner data. Int J Remote Sens. 1985;6(1):115-125.

mote Sens. 2009;30(21):5737-5750.

71-82.

196 Environmental Applications of Remote Sensing

Oceans are a vast, complex world where underwater sound is the most efficient tool available to understand its detailed characteristics. However the underwater chan‐ nel has a very complex geometrical and material structure and hence special techni‐ ques are required to model it. Analytical solutions are feasible only when one makes gross assumptions and approximations. Several numerical and semi-numeri‐ cal 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 ex‐ plore in this chapter the scope and feasibility of finite element method in underwa‐ ter 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 ei‐ gensolution for this channel is required for imposing the radiation condition at the truncation boundary. We have cast this eigenproblem in a variational form and em‐ ployed a Rayleigh-Ritz method to obtain an approximate eigensolution. This ap‐ proach has provided a good approximation of the depth eigenmodes in a compact semi-analytic form. We have employed our finite element algorithm to model sever‐ al range- and depth-dependent ocean problems. Our numerical study has establish‐ ed that our finite element algorithm gives accurate results with reasonable effort. In particular, our finite element approach is most appropriate for shallow water prob‐ lems 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

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

to several other ocean acoustic problems encountered in the remote sensing of ocean environment.

**Keywords:** Wave propagation, scattering, ocean wave guide, irregular boundaries

### **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 sound1 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 transmis‐ sion, 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 shad‐ owing 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

<sup>1</sup> 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.

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**

to several other ocean acoustic problems encountered in the remote sensing of

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 sound1 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

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 transmis‐ sion, 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

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 shad‐ owing 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

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

1 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

**Keywords:** Wave propagation, scattering, ocean wave guide, irregular boundaries

[1-5] has been driven by applications more than other branches of ocean science.

significant acoustic fluctuations over long transmission paths.

ocean environment.

198 Environmental Applications of Remote Sensing

**1. Introduction**

channel.

underwater channels.

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 rangedependent and range-independent problems, but are most commonly used for rangeindependent 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 hori‐ zontal. 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 com‐ putational 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 computa‐ tionally 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**

This approximation turns out to include negligible errors beyond a wavelength or so from the

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 hori‐ zontal. 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

source.

200 Environmental Applications of Remote Sensing

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, heteroge‐ neous 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**

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, heteroge‐ neous 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

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

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

waveguide.

202 Environmental Applications of Remote Sensing

form.

sensing of ocean environment.

The fluid domain Ω=Ω*<sup>I</sup>* + Ω*O* (Fig. 1) of the waveguide problem consists of the inner domain Ω*I* truncated by the artificial radiation boundary *SR*, and the outer domain Ω*O* (far-field domain). The waveguide is assumed to be axially symmetric about the vertical axis containing a source at depth *zs*, 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 (*SF* ), and at the bottom by a seabed of arbitrary topography (*SB*). The waveguide is assumed to have unbounded range. For timeharmonic linear acoustic waves with the pressure field denoted as *p* ^(*<sup>r</sup>*, *<sup>z</sup>*, *<sup>t</sup>*)= *<sup>p</sup>*(*r*, *<sup>z</sup>*)*<sup>e</sup>* <sup>−</sup>*iω<sup>t</sup>* , *ω* 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\_\circ \delta(r) \delta\left(z - z\_s\right),
\tag{1}
$$

where ∇ is the gradient operator, *ρ* the density of the acoustic fluid, *k* the acoustic wavenum‐ ber, *c* the local speed of sound, and *f <sup>o</sup>* defines the point source at *r* =0 and *z* = *zs*.

**Figure 1.** Geometry of the ocean waveguide

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

$$p = 0 \qquad \text{on } S\_{\mathbb{F}} \tag{2}$$

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 \qquad \text{on } S\_{\mathbb{B}'} \tag{3}$$

where *S*B denotes the sea bottom.

For the purpose of FE modeling, the waveguide, which is unbounded in range, is truncated at *r* =*rb*, and the truncation boundary is treated as the radiation boundary *SR*, 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 \mathbf{n}} + \alpha p\_m = 0, \qquad m = 1, 2, \dots, M \quad \text{on } S\_{\mathbb{R}^\nu} \tag{4}$$

where *M* denotes the number of propagating modes, and the damper coefficient *αm* associated with the *m-*th mode is given by

1 <sup>2</sup> , *<sup>m</sup> <sup>r</sup> rm* a= - *ik* (5)

where *krm* 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 *SR*, acoustic pressure may be expressed as a sum of normal pressure modes as

1 ( ) ( ), *M m m pz p z* = <sup>=</sup> å (6)

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} + \sum\_{m=1}^{M} \alpha\_m p\_m \quad \text{on } \mathbb{S}\_{\mathbb{R}}\,. \tag{7}$$

Denoting a normal-mode function at the radiation boundary by *f <sup>m</sup>*(*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), \qquad m = 1, 2, \cdots, M \,\,\,\tag{8}$$

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

$$p(\mathbf{z}) + \sum\_{m=1}^{M} a\_m \alpha\_m f\_m(\mathbf{z}) = 0 \qquad \text{on } S\_{\mathbb{R}^\times} \tag{9}$$

where the constants *am* are determined by using the (1 / *ρ*(*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*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 *am*.

#### **5. Constraints**

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

$$\left[\mathbf{C}\right]\left(p(z), a\_1, a\_2, \dots, a\_M\right)^T = 0 \quad \text{on } S\_{\mathbb{R}'} \tag{10}$$

where

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

For the purpose of FE modeling, the waveguide, which is unbounded in range, is truncated at *r* =*rb*, and the truncation boundary is treated as the radiation boundary *SR*, on which a suitable approximate radiation condition should be imposed. Here, the boundary damper approach

¶ <sup>=</sup> ¶ (3)

¶ <sup>L</sup> (4)

= - *ik* (5)

<sup>=</sup> å (6)

¶ å (7)

( ) ( ), 1,2, , , *m mm p z af z m M* = = L (8)

<sup>B</sup> 0 on , *<sup>p</sup> <sup>S</sup>*

[26] has been adopted. The first-order cylindrical damper equation may be written as

0, 1,2, , on , *<sup>m</sup> m R <sup>p</sup> p m MS*

> 1 <sup>2</sup> , *<sup>m</sup> <sup>r</sup> rm*

where *M* denotes the number of propagating modes, and the damper coefficient *αm* associated

where *krm* 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 *SR*, acoustic

> 1 ( ) ( ), *M m m pz p z* =

> > 1

=

a

*<sup>p</sup> p S*

*M*

*m*

*r*

¶ +

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

*mm R*

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

on .

*n*

+= =

a

pressure may be expressed as a sum of normal pressure modes as

boundary. In other words,

204 Environmental Applications of Remote Sensing

where *S*B denotes the sea bottom.

with the *m-*th mode is given by

written, using Eqs. (4) and (6), as

*n* a

¶

$$\mathbf{f}\left[\mathbf{C}\right] = \left[\mathbf{1}\_{\prime} - f\_1(\mathbf{z})\_{\prime} - f\_2(\mathbf{z})\_{\prime}, \dots, -f\_M(\mathbf{z})\_{\prime}\right],\tag{11}$$

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 Ω*<sup>S</sup>* surrounding the source has been excluded so that the computational domain is Ω¯ *<sup>I</sup>* =Ω*<sup>I</sup>* −Ω*<sup>S</sup>* . Consider the following axisymmetric functional *I*(*p*) defined in the cylindrical coordinate system (*r*, *z*) (see Fig. 1):

$$I(p) = \frac{1}{2} \int\_{\overline{\Omega}\_l} \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 \, dr \, d\overline{z} + \frac{1}{2} \sum\_{m=1}^M \int\_{S\_{\mathbb{R}}} \frac{1}{\rho} \alpha\_m p\_m^2 r d\overline{z} - \int\_{S\_{\mathbb{D}} \times \mathbb{S}\_N} \frac{1}{\rho} \frac{\partial p}{\partial n} p d\mathcal{S},\tag{12}$$

where *S*D denotes the surface on which a Dirichlet boundary condition is prescribed and *S*<sup>N</sup> the surface with prescribed Neumann boundary condition, and the other domains of integra‐ tion are identified in Fig. 1.

It can readily be shown that the variational condition

$$
\delta I = 0 \tag{13}
$$

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 Ω¯ *<sup>I</sup>* (which excludes the source) of the axisymmetric waveguide in Fig. 1 may be discretized using eight-noded axisymmetric quadrilateral elements with *Co* 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 \sum\_{j=1}^{\bar{n}} \overline{p}\_{e\_j} N\_{\backslash}(\xi, \zeta) = \left[N\right]^T \left\{\overline{p}\_e\right\},\tag{14}$$

where *n*˜ denotes the number of element nodes (eight in the present study), *p*¯*ej* the nodal pressure variable/degrees of freedom (dofs) and *Nj* (*ξ*, *ζ*) the polynomial shape function in the parametric coordinates (*ξ*, *ζ*) 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(\boldsymbol{p}\_{\boldsymbol{\epsilon}}) \approx \frac{1}{2} \left\{ \overline{\boldsymbol{p}}\_{\boldsymbol{\epsilon}} \right\}^{T} \left( \left[ \boldsymbol{K}\_{\boldsymbol{\epsilon}} \right] - \left[ \boldsymbol{M}\_{\boldsymbol{\epsilon}} \right] \right) \left\{ \overline{\boldsymbol{p}}\_{\boldsymbol{\epsilon}} \right\} + \frac{1}{2} \sum\_{m=1}^{M} \left\{ \overline{\boldsymbol{p}}\_{cm} \right\}^{T} \left[ \overline{\boldsymbol{R}}\_{cm} \right] \left\{ \overline{\boldsymbol{p}}\_{cm} \right\} - \left\{ \overline{\boldsymbol{p}}\_{\boldsymbol{\epsilon}} \right\}^{T} \left\{ f\_{\boldsymbol{\epsilon}} \right\},\tag{15}$$

where {*p*¯*em*} 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*(*pe*) 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* ′ may be defined as

$$I'=I+\frac{1}{2}\left\{\overline{p}\_e'\right\}^T \left[\overline{\mathbb{C}\_e}\right]^T \left[\overline{\mathcal{J}}\_P\right] \left[\overline{\mathbb{C}\_e}\right] \left\{\overline{p}\_e'\right\}\_{'} \tag{16}$$

where *Ce* denotes the constraint matrix in Eq. (11) specific to an element. The penalty coefficient matrix *βP* above may be chosen to be diagonal for convenience, with *βPm* denoting the penalty parameter associated with the *m-*th mode. Equation (16) may be expanded as

$$I'(p'\_{\varepsilon}) = \frac{1}{2} \left(\overline{p}'\_{\varepsilon}\right)^{T} \left(\left[\overline{\mathbf{K}}'\_{\varepsilon}\right] - \left[\overline{\mathbf{M}}'\_{\varepsilon}\right] + \left[\overline{\mathbf{R}}'\_{\varepsilon}\right] + \left[\mathbf{C}\_{\varepsilon}\right]^{T} \left[\left[\mathcal{B}\_{p}\right]\right] \overline{\mathbf{C}}\_{\varepsilon}\right) \left\{\overline{p}'\_{\varepsilon}\right\} - \left\{\overline{p}'\_{\varepsilon}\right\}^{T} \left\{f'\_{\varepsilon}\right\} \overline{\mathbf{C}}\_{\varepsilon}$$

where the enlarged element dof vector is defined as

( ) ( )

It can readily be shown that the variational condition

*I*

206 Environmental Applications of Remote Sensing

tion are identified in Fig. 1.

**7. Finite element model**

variable *p* may then be written as

The finite fluid domain Ω¯

r

2 2

R

ì ü = +- + í ý - î þ <sup>ò</sup> åò ò (12)

r

a

1

*<sup>M</sup> p p <sup>p</sup> r z m m <sup>n</sup> <sup>S</sup> <sup>m</sup> S S*

 ¶ ¶ ¶ ¶ ¶ ¶

where *S*D denotes the surface on which a Dirichlet boundary condition is prescribed and *S*<sup>N</sup> the surface with prescribed Neumann boundary condition, and the other domains of integra‐

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).

1 may be discretized using eight-noded axisymmetric quadrilateral elements with *Co* 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

> (,) ( , ) , *<sup>n</sup> <sup>T</sup> ej j e*

where *n*˜ denotes the number of element nodes (eight in the present study), *p*¯*ej* the nodal

parametric coordinates (*ξ*, *ζ*) 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:

{ } ( ){ } { } { } { }{ } 1 1

=

1 ( ) , *<sup>M</sup> <sup>T</sup> T T e e e ee em em em e e m Ip p K M p p R p p f*

where {*p*¯*em*} 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*(*pe*) 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

» -+ - é ùé ù é ù ë ûë û å ë û (15)

» = é ù å ë û

*prz p N N p* x z

1

*j*

pressure variable/degrees of freedom (dofs) and *Nj*

2 2

=

%

2 2 1 1 2 2 1 1 <sup>2</sup> <sup>1</sup>

*I p k p r dr dz p rdz pdS*

W = +

d

( ) ,

D N

*I* = 0 (13)

*<sup>I</sup>* (which excludes the source) of the axisymmetric waveguide in Fig.

{ }

(*ξ*, *ζ*) the polynomial shape function in the

(14)

 r

$$\left\{ \left\{ \overline{p}'\_e \right\} = \begin{Bmatrix} \left\{ \overline{p}\_e \right\} \\ \left\{ a \right\} \end{Bmatrix} \right\}. \tag{17}$$

(19)

The enlarged stiffness, mass and damping matrices, and load vector in expanded Eq. (16), consistent with {*p*¯′ *e* } in Eq. (17), are given by

$$
\begin{bmatrix} \mathbf{K}'\_{\epsilon} \\ \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} \mathbf{K}\_{\epsilon} \\ 0 & 0 \end{bmatrix} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \qquad \begin{bmatrix} \mathbf{M}'\_{\epsilon} \\ \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} \mathbf{M}\_{\epsilon} \\ 0 & 0 \end{bmatrix} & \begin{bmatrix} \mathbf{R}'\_{\epsilon} \\ \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \begin{bmatrix} \mathbf{R}\_{\epsilon} \\ \end{bmatrix} \end{bmatrix} \qquad \begin{Bmatrix} f'\_{\epsilon} \end{Bmatrix} = \begin{Bmatrix} \begin{Bmatrix} \langle f\_{\epsilon} \end{Bmatrix} \\ \mathbf{0} \end{Bmatrix}. \tag{18}
$$

The matrices *Ke* , *Me* , and *Re* in Eq. (18) are traditionally called the element stiffness, mass and radiation damping matrices, and { *f <sup>e</sup>*} the load vector, respectively. They are given as follows:

$$
\left[\begin{array}{c} \mathbf{K}\_c \\ \end{array}\right] = \int\_{\Omega\_c} \frac{1}{\rho} \left[\begin{array}{c} \nabla \mathbf{N} \\ \end{array}\right]^T \left[\begin{array}{c} \nabla \mathbf{N} \\ \end{array}\right] d\Omega \tag{\text{a}}
$$

$$
\left[\begin{array}{c} \boldsymbol{M}\_{\boldsymbol{\epsilon}} \end{array}\right] = \int\_{\boldsymbol{\Omega}\_{\boldsymbol{\epsilon}}} \frac{k^{2}}{\rho} \left[\boldsymbol{N}\right]^{T} \left[\boldsymbol{N}\right] d\boldsymbol{\Omega} \tag{b}
$$

$$\left[\boldsymbol{\mathcal{R}}\_{\boldsymbol{\epsilon}}\right] = \text{diag}\left[\boldsymbol{\mathcal{R}}\_{\boldsymbol{\epsilon}1'}\boldsymbol{\mathcal{R}}\_{\boldsymbol{\epsilon}2'}, \dots, \boldsymbol{\mathcal{R}}\_{\boldsymbol{\epsilon}\boldsymbol{\mathcal{M}}}\right] \tag{3}$$

$$\mathcal{R}\_{cm} = \left\{ f\_{zw} \right\}^T \left[ \overline{\mathcal{R}}\_{cm} \right] \left\{ f\_{zw} \right\} \qquad \left[ \overline{\mathcal{R}}\_{cm} \right] = \mathcal{R}\_{m} \int\_{S\_{R}} \frac{1}{\rho} \left[ \boldsymbol{N} \right]^T \left[ \boldsymbol{N} \right] dS \quad \text{(d)}$$

$$\left\{ f\_{zw} \right\} = \left\{ f\_{zw}(\mathbf{z}\_1)\_{\prime} f\_{zw}(\mathbf{z}\_2)\_{\prime} \dots f\_{zw}(\mathbf{z}\_n) \right\}^T \tag{6}$$
 
$$\mathbf{f}\_{zw} \dots \mathbf{f}\_{zw} \mathbf{z}\_T \dots \mathbf{z}\_T \dots$$

$$\left\{f\_{\epsilon}\right\} = \int\_{\mathcal{S}\_{\mathsf{M}}} \frac{1}{\rho} p\_v \left[\!\!\!/\!\!N\right]^\top dS \tag{f}$$

where *N* denotes the shape-function matrix [see Eq. (14)], *p<sup>ν</sup>* =∂ *p* / ∂*n*, and *f zm*(*zj* ) denotes the *j*-th nodal value of the *m-*th mode on a finite element in contact with the radiation boundary *SRe*. The steps required to derive Eq. (19c) are outlined in Appendix A. It is of interest to note that the radiation-damping matrix *Re* in Eq. (19c) implies uncoupled modal participation. However, the constraint term involving the matrix *Ce* 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* ′ in expanded Eq. (16) is obtained by setting

$$\frac{\partial I'}{\partial \left\{ \overline{p}'\_{\epsilon} \right\}} = 0 \tag{20}$$

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

$$
\left( \left[ \boldsymbol{K}'\_{\boldsymbol{\epsilon}} \right] - \left[ \boldsymbol{M}'\_{\boldsymbol{\epsilon}} \right] + \left[ \boldsymbol{R}'\_{\boldsymbol{\epsilon}} \right] + \left[ \boldsymbol{C}\_{\boldsymbol{\epsilon}} \right]^{T} \left[ \boldsymbol{\mathcal{J}} \boldsymbol{\mathcal{J}}\_{\mathbb{P}} \right] \left[ \boldsymbol{C}\_{\boldsymbol{\epsilon}} \right] \right) \left\{ \overline{p}'\_{\boldsymbol{\epsilon}} \right\} = \left\{ f'\_{\boldsymbol{\epsilon}} \right\}.\tag{21}
$$

It may be noted that if the penalty matrix *β<sup>P</sup>* =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 *R* ′ *<sup>e</sup>* and the constraint matrix *Ce* 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[ \boldsymbol{K}\_{e} \right] - \left[ \boldsymbol{M}\_{e} \right] \right) \left( \overline{p}\_{e} \right) = \left\{ f\_{e} \right\}. \tag{22}$$

The radiation-damping matrix *R* ′ *<sup>e</sup>* 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 *SR*. Carrying out thus the finite element assemblage operation, yields the following global finite element equations:

$$\left( \left[ \begin{array}{c} \text{K}' \end{array} \right] - \left[ \begin{array}{c} \text{M}' \end{array} \right] + \left[ \begin{array}{c} \text{R}' \end{array} \right] + \left[ \begin{array}{c} \text{C}' \end{array} \right] \right) \left\{ \overline{p}' \right\} = \left\{ f' \right\}, \tag{23}$$

where the global solution vector {*p*¯′ } consists of all the pressure dof in the computational domain as well as the unknown vector {*a*} in Eq. (8).

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

$$
\begin{aligned}
\left[\begin{array}{c}\mathbf{K}'\right] = \sum\limits\_{\epsilon} \begin{bmatrix} K'\_{\epsilon} \end{bmatrix} \qquad \left[\begin{array}{c}\mathbf{M}' \end{array}\right] = \sum\limits\_{\epsilon} \begin{bmatrix} M'\_{\epsilon} \end{bmatrix} \qquad \left[\begin{array}{c}\mathbf{R}' \end{array}\right] = \sum\limits\_{\epsilon} \begin{bmatrix} R'\_{\epsilon} \end{bmatrix} \\\ \left[\begin{array}{c}\mathbf{C}' \end{array}\right] = \sum\limits\_{\epsilon} \begin{bmatrix} \mathbf{C}\_{\epsilon} \end{bmatrix}^{T} \begin{bmatrix} \mathcal{J}\_{p} \end{bmatrix} \begin{bmatrix} \mathbf{C}\_{\epsilon} \end{bmatrix} \qquad \left\{\begin{array}{c}f' \end{array}\right\} = \sum\limits\_{\epsilon} \begin{Bmatrix} f'\_{\epsilon} \end{Bmatrix}\tag{24}
$$

where Σ*e* denotes the standard finite element assemblage operation [15].

#### **8. Modeling of a point source**

where *N* denotes the shape-function matrix [see Eq. (14)], *p<sup>ν</sup>* =∂ *p* / ∂*n*, and *f zm*(*zj*

The stationary condition of the potential *I* ′

208 Environmental Applications of Remote Sensing

The radiation-damping matrix *R* ′

where the global solution vector {*p*¯′

domain as well as the unknown vector {*a*} in Eq. (8).

element equations:

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

*j*-th nodal value of the *m-*th mode on a finite element in contact with the radiation boundary *SRe*. The steps required to derive Eq. (19c) are outlined in Appendix A. It is of interest to note that the radiation-damping matrix *Re* in Eq. (19c) implies uncoupled modal participation. However, the constraint term involving the matrix *Ce* in expanded Eq. (16) brings about modal coupling. The various integrals above are defined over relevant finite element domains.

> { } <sup>0</sup> *e I p*

( ){} {}. *<sup>T</sup> K M R C Cp f e e e e Pe e e* é ù é ù é ù é ù é ùé ù ¢ ¢¢ - ++

choice of the penalty parameters is essential. The radiation-damping matrix *R* ′

for this case, the FE equation may be deduced from Eqs. (18) and (21) as

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

b

It may be noted that if the penalty matrix *β<sup>P</sup>* =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

constraint matrix *Ce* in Eq. (21) correspond to elements on the radiation boundary. Hence,

only for elements that share one or more of their boundaries with the artificial boundary *SR*. Carrying out thus the finite element assemblage operation, yields the following global finite

> *e ee e ee T*

*C C Cf f* b

é ù é ù é ùé ù ¢ = = ¢ ¢ ë û ë û ë ûë û

*K K M MR R*

éù é ù é ù é ù éù é ù ¢ ¢ ¢ ¢¢ ¢ = == ëû ë û ë û ë û ëû ë û

å åå

*e Pe e e e*

¢ ¢ = ë û ë û ë û ë û ë ûë û (21)

( ){} {}. *K Mp f e ee e* é ùé ù - = ë ûë û (22)

(éùé ùéùéù *K M R Cp f* ¢ ¢ ¢ ¢¢ ¢ - ++ = ){} {}, ëûë ûëûëû (23)

{} {}

å å (24)

*<sup>e</sup>* in Eq. (19c), which is complex in view of Eq. (5), is defined

} consists of all the pressure dof in the computational

in expanded Eq. (16) is obtained by setting

¶ ¢ <sup>=</sup> ¶ ¢ (20)

) denotes the

*<sup>e</sup>* and the

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

$$\begin{bmatrix} A \\ \end{bmatrix} \begin{Bmatrix} \overline{p}' \end{Bmatrix} = \begin{Bmatrix} f' \end{Bmatrix} \tag{25}$$

It may be noted that for an acoustic medium with real sound speed, the coefficient matrix *A* above is complex even though *K* ′ , *M* ′ , and *C* ′ are real. This is because { *f* ′ } is complex. Also, note that *<sup>R</sup>* ′ is complex due to the presence of *αm* in Eq. (19d). For a lossy medium modeled with complex sound speed, *M* ′ is also complex. Although *A* 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 *α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 *A* in Eq. (25); i.e.,

$$\beta\_{\text{Pu}} = \frac{\beta\_s}{n'} \sum\_{i=1}^{n'} \left| A\_{ii} \right| \tag{26}$$

where *n* ′ denotes the total number of FE equations/dof and *βs* a user-specified penalty scale factor. Computations indicate that the results are stable over a wide range of *β<sup>s</sup>* values. The results reported here have been obtained using *β<sup>s</sup>* =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 rangedependent. 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 rangedependent 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**

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

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 *α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

1

factor. Computations indicate that the results are stable over a wide range of *β<sup>s</sup>* values. The

The sound speed in an ocean-acoustic waveguide is in general both depth- and rangedependent. 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 rangedependent 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

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

*n s Pm ii i A*

¢ =

*n* b

**10. Normal modes in an ocean waveguide with depth dependence**

b

,

denotes the total number of FE equations/dof and *βs* a user-specified penalty scale

<sup>=</sup> ¢ å (26)

they hold promise for high frequency FE models.

results reported here have been obtained using *β<sup>s</sup>* =100.

this aspect with specific reference to shallow-water waveguides.

coefficient matrix *A* in Eq. (25); i.e.,

210 Environmental Applications of Remote Sensing

where *n* ′

For the cylindrically symmetric waveguide having depth-dependent density *ρ* 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{\alpha^2}{c^2(z)}p = -\frac{\delta(r)\delta(z-z\_s)}{2\pi r},\tag{27}$$

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* = *zs*, with *δ* 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) \tag{28}$$

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{\mathcal{R}}(r)}{dr}\right) + k\_r^2 \overline{\mathcal{R}}(r) = 0\tag{29}$$

$$
\rho \rho(z) \frac{d}{dz} \left( \frac{1}{\rho(z)} \frac{dZ(z)}{dz} \right) + \left( \frac{\alpha^2}{c^2(z)} - k\_r^2 \right) Z(z) = 0 \tag{30}
$$

where *ω* denotes the circular frequency and *kr* 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 *R*¯(*r*), and Eq. (30) pertains to the depth modes *Z*(*z*). Choosing a pressure release boundary at the top (*z* =0) and a mixed/Robin boundary condition at the seabed (*z* =*D*1), the boundary conditions of our problem are written as [1, 30, 33]

$$Z(0) = 0\tag{31}$$

$$Z(D\_1) + \frac{g\left(k\_r^2\right)}{\rho} \frac{dZ(D\_1)}{dz} = 0\tag{32}$$

where

$$\lg(k\_r^2) = \rho\_b / \sqrt{\left(k\_r^2 - \frac{\alpha^2}{c^2}\right)}$$

with *ρ<sup>b</sup>* 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 *kr* 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 halfspace below the water column are also considered here [1]. The governing equation for the waves in this fluid half-space is given as,

$$\rho\_{\rm b}(\mathbf{z}) \frac{d}{dz} \left( \frac{1}{\rho\_{\rm b}(\mathbf{z})} \frac{dZ\_{\rm b}(\mathbf{z})}{dz} \right) + \left( \frac{\alpha^2}{c\_{\rm b}^2(\mathbf{z})} - k\_r^2 \right) Z\_{\rm b}(\mathbf{z}) = 0, \qquad D\_1 \le \mathbf{z} \le \infty \tag{33}$$

where *Z*<sup>b</sup> denotes the depth function in the fluid half-space having depth-dependent density *ρ*<sup>b</sup> and sound speed *c*b. The interface conditions at the seabed are given by the kinematic and dynamic conditions,

$$\begin{aligned} Z(D\_1) &= Z\_b(D\_1) & \text{(a)}\\ \frac{1}{\rho} \frac{dZ(D\_1)}{dz} &= \frac{1}{\rho\_b} \frac{dZ\_b(D\_1)}{dz} & \text{(b)} \end{aligned} \tag{34}$$

In addition, the depth mode *Z*<sup>b</sup> should remain bounded as *z* →*∞*. 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 *Π*(*Z*) and *Πb*(*Zb*) defined respectively in the water column and the half-space as

$$\prod(Z) = \frac{1}{2} \int\_0^{D\_1} \left[ \frac{1}{\rho(z)} \left( \frac{dZ}{dz} \right)^2 - \frac{\alpha^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 \Big|\_0^{D\_1} \tag{35}$$

$$\begin{split} \Pi\_{b}(\boldsymbol{Z}\_{b}) &= \frac{1}{2} \int\_{D\_{1}}^{D\_{2}} \left[ \frac{1}{\rho\_{b}(\boldsymbol{z})} \left( \frac{d\boldsymbol{Z}\_{b}}{d\boldsymbol{z}} \right)^{2} - \frac{\alpha^{2}}{\rho\_{b}(\boldsymbol{z})c^{2}(\boldsymbol{z})} \boldsymbol{Z}\_{b}^{2} + \frac{1}{\rho\_{b}(\boldsymbol{z})} \boldsymbol{k}\_{r}^{2} \boldsymbol{Z}\_{b}^{2} \right] d\boldsymbol{z} \\ &- \frac{1}{\rho\_{b}(\boldsymbol{z})} \boldsymbol{Z}\_{b\upsilon}(\boldsymbol{D}\_{1}) \boldsymbol{Z}\_{b}(\boldsymbol{D}\_{1}) + \frac{1}{2} \beta \frac{\boldsymbol{Z}\_{b}^{2}(\boldsymbol{D}\_{2})}{\rho\_{b}}, \qquad \boldsymbol{D}\_{2} \to \infty \end{split} \tag{36}$$

where suffix *ν* 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\_{\boldsymbol{\nu}}(D\_1) = \frac{1}{\rho\_b} \frac{dZ\_{\boldsymbol{\nu}}(D\_1)}{dz} \,, \tag{37}$$

and in Eq. (36),

1 () <sup>2</sup> () 0 *<sup>r</sup> d dR r r kRr*

*z k Zz dz z dz c z*

( ) <sup>2</sup>

( ) ( ) <sup>0</sup> *<sup>r</sup> g k dZ D Z D* r*dz*

> 2 2 <sup>2</sup> () / *rb r gk k*

with *ρ<sup>b</sup>* 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-

[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 halfspace below the water column are also considered here [1]. The governing equation for the

æ ö æ ö ç ÷ + - = £ £¥ ç ÷

where *Z*<sup>b</sup> denotes the depth function in the fluid half-space having depth-dependent density *ρ*<sup>b</sup> and sound speed *c*b. The interface conditions at the seabed are given by the kinematic and

2 b 2 b 2 b 1

w

*z k Zz D z*

æ ö = - ç ÷

r

1

2

*c* w

è ø

ç ÷ + =

2

2 1 () ( ) () 0 ( ) ( ) *<sup>r</sup>*

æ ö æ ö ç ÷ +- = ç ÷ è ø è ø

2

the square of the radial/horizontal wavenumber. Eq. (29) evidently pertains to the radial/ horizontal modes *R*¯(*r*), and Eq. (30) pertains to the depth modes *Z*(*z*). Choosing a pressure release boundary at the top (*z* =0) and a mixed/Robin boundary condition at the seabed

w

2

è ø (29)

, the separation constant, which turns out to be

*Z*(0) 0 = (31)

+ = (32)

2

. Porter and Reiss

(33)

(30)

æ ö

*r dr dr*

(*z* =*D*1), the boundary conditions of our problem are written as [1, 30, 33]

1

Liouville problem because Eq. (32) depends on the unknown eigenvalue *kr*

<sup>1</sup> ( ) ( ) ( ) 0, ( ) ( ) *<sup>r</sup>*

è ø è ø

b b

*dz z dz c z*

*d dZ z*

r

*d dZ z*

r

r

212 Environmental Applications of Remote Sensing

where *ω* denotes the circular frequency and *kr*

waves in this fluid half-space is given as,

r

dynamic conditions,

where

$$\frac{1}{\rho\_b} Z\_{\rm by}(D\_1) = \frac{1}{\rho} \frac{dZ(D\_1)}{dz} \tag{38}$$

In addition, Eq. (34a) should be imposed. Then, it can be easily shown that the variational condition *δΠ* =0 leads to Eq. (30) and the boundary conditions in Eq. (31) as well as the interface condition in Eq. (37), where *δ* denotes the first variation. Similarly, the variational condition *δΠ<sup>b</sup>* =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{dZ\_{\flat}(D\_{2})}{dz} + \beta Z\_{\flat}(D\_{2}) = 0,\tag{39}$$

where (see Eq. (32)),

$$\beta = \sqrt{\left(k\_r^2 - \frac{\alpha^2}{c\_b^2}\right)}$$

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

Case 1: *D*2 is finite, *β* →0

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

Case 2: *β* is finite, *D*2→*∞*

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

Case 3: *D*2 and *β* 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{aligned} Z(\mathbf{z}) &\approx \sum\_{j=1}^{n\_1} \overline{\phi}\_j \mathbb{1}\_{\boldsymbol{\uprho}}(\mathbf{z}) \qquad 0 \le \mathbf{z} \le D\_1 \quad \text{(a)}\\ Z\_{\boldsymbol{\uprho}}(\mathbf{z}) &\approx \sum\_{j=1}^{n\_2} \overline{\phi}\_{\boldsymbol{\uprho}} \mathbb{1}\_{\boldsymbol{\uprho}}(\mathbf{z}) \quad D\_1 \le \mathbf{z} \le \boldsymbol{\uprho} \quad \text{(b)}\end{aligned} \tag{40}$$

where *n* =*n*<sup>1</sup> + *n*2, and *ψ<sup>j</sup>* and *ψbj* denote the known mode function (coordinate function) satisfying the kinematic boundary condition in the water column and *ϕ*¯ *<sup>j</sup>* 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(\mathbf{z}) \approx \sum\_{j=1}^{n} \overline{\phi}\_{j} \mathbb{1}\_{j}(\mathbf{z}) \quad 0 \le \mathbf{z} \le D\_{\mathbf{z}} \tag{41}$$

where we have combined the depth modes of a two-layer isovelocity waveguide as one combined set with redefined coefficients *ϕ*¯. 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{\Pi} = \Pi + \Pi\_b = \frac{1}{2} \left\{ \overline{f} \right\}^T \left[ \mathbf{K}^\prime \right] \left\{ \overline{f} \right\} - \frac{1}{2} \alpha^2 \left\{ \overline{f} \right\}^T \left[ \mathbf{K}^{\text{II}} \right] \left\{ \overline{f} \right\} + \frac{1}{2} k\_r^2 \left\{ \overline{f} \right\}^T \left[ \mathbf{K}^{\text{III}} \right] \left\{ \overline{f} \right\} \tag{42}$$

where

In addition, Eq. (34a) should be imposed. Then, it can be easily shown that the variational condition *δΠ* =0 leads to Eq. (30) and the boundary conditions in Eq. (31) as well as the interface condition in Eq. (37), where *δ* denotes the first variation. Similarly, the variational condition *δΠ<sup>b</sup>* =0 leads to Eq. (33), and the interface conditions in Eq. (34a) and Eq. (38). In addition, at

> 2 ( ) ( ) 0, *<sup>b</sup> b*

> > 2

*c* w

(39)

+ = b

> 2 *r* 2 *b*

æ ö = - ç ÷ è ø

*k*

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

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

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

1

( ) ( ) 0 (a)

() () (b)

» £ £¥

1

and *ψbj* denote the known mode function (coordinate function)

(40)

an unknown constant,

2

*dz*

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

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

seabed of infinite extent with uniform sound speed and density.

1

*n*

*j n b bj bj j*

=

å

2

1

satisfying the kinematic boundary condition in the water column and *ϕ*¯ *<sup>j</sup>*

=

å

1

*j*

f y

f y

*j*

*Zz z D z*

*Zz z z D*

» £ £

*dZ D Z D*

b

*z* =*D*2, we obtain the condition

214 Environmental Applications of Remote Sensing

where (see Eq. (32)),

Case 1: *D*2 is finite, *β* →0

Case 2: *β* is finite, *D*2→*∞*

Case 3: *D*2 and *β* are finite

where *n* =*n*<sup>1</sup> + *n*2, and *ψ<sup>j</sup>*

form

by a rigid boundary.

$$\begin{aligned} K\_{\boldsymbol{\uprho}}^{\boldsymbol{D}} &= \int\_{0}^{\boldsymbol{D}\_{1}} \frac{1}{\rho(\boldsymbol{z})} \frac{d\boldsymbol{\uprho}\_{1}}{dz} \frac{d\boldsymbol{\uprho}\_{1}}{dz} dz + \int\_{\boldsymbol{D}\_{1}}^{\boldsymbol{D}\_{2}} \frac{1}{\rho\_{\boldsymbol{b}}(\boldsymbol{z})} \frac{d\boldsymbol{\uprho}\_{\boldsymbol{b}}}{dz} \frac{d\boldsymbol{\uprho}\_{\boldsymbol{b}}}{dz} dz \\\ K\_{\boldsymbol{\uprho}}^{\boldsymbol{D}} &= \int\_{0}^{\boldsymbol{D}\_{1}} \frac{1}{\rho(\boldsymbol{z})c^{2}(\boldsymbol{z})} \boldsymbol{\uprho}\_{\boldsymbol{\uprho}} \boldsymbol{\uprho}\_{\boldsymbol{b}} dz + \int\_{\boldsymbol{D}\_{1}}^{\boldsymbol{D}\_{2}} \frac{1}{\rho\_{\boldsymbol{b}}(\boldsymbol{z})c\_{\boldsymbol{b}}^{2}(\boldsymbol{z})} \boldsymbol{\uprho}\_{\boldsymbol{b}} \boldsymbol{\uprho}\_{\boldsymbol{b}} dz \\\ K\_{\boldsymbol{\uprho}}^{\boldsymbol{D}\_{1}} &= \int\_{0}^{\boldsymbol{D}\_{1}} \frac{1}{\rho(\boldsymbol{z})} \boldsymbol{\uprho}\_{\boldsymbol{b}} \boldsymbol{\uprho}\_{\boldsymbol{b}} dz + \int\_{\boldsymbol{D}\_{1}}^{\boldsymbol{D}\_{2}} \frac{1}{\rho\_{\boldsymbol{b}}(\boldsymbol{z})} \boldsymbol{\uprho}\_{\boldsymbol{b}} \boldsymbol{\uprho}\_{\boldsymbol{b}} dz \end{aligned} \tag{43}$$

It has been assumed in the above that the contribution due to the term in Eq. (36) is negligible as *D*2→*∞*. 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 \phi\_{\rangle}} = 0, \quad j = 1, 2, \dots, n \tag{44}$$

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

$$\mathbb{E}\left(\boldsymbol{\phi}^{2}\left[\boldsymbol{\K}^{\mathrm{II}}\right] - \left[\boldsymbol{\K}^{\mathrm{I}}\right]\right)\left\{\overline{\boldsymbol{\phi}}\right\} = k\_{r}^{2}\left[\boldsymbol{\K}^{\mathrm{III}}\right]\left\{\overline{\boldsymbol{\phi}}\right\}\tag{45}$$

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

$$\left(k\_{\gamma'}^2 \left\{\phi^{(\prime)}\right\}\right), \qquad \qquad j = 1, 2, \dots \tag{46}$$

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 *krj* 2 and the eigenvectors {*ϕ*¯( *<sup>j</sup>*) }, the eigenfunctions /depth modes may be written, using Eq. (40), as

$$\text{var}\_{\rangle}(z) = \left\{ \overline{\Phi}^{(\prime)} \right\}^{T} \left\{ \psi(z) \right\} \tag{47}$$

where

$$\left\{\boldsymbol{\psi}(\boldsymbol{z})\right\}^{T} = \left\{\boldsymbol{\psi}\_{1}(\boldsymbol{z}), \boldsymbol{\psi}\_{2}(\boldsymbol{z}), \dots, \boldsymbol{\psi}\_{n}(\boldsymbol{z})\right\} \tag{48}$$

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 approx‐ imating 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 KIII -orthogonal, it can easily be shown that the eigenfunctions in Eq. (47) satisfy the following orthogonality condition:

$$\int\_{0}^{D\_{1}} \frac{1}{\rho(z)} Z\_{i} Z\_{j} dz + \int\_{D\_{1}}^{D\_{2}} \frac{1}{\rho\_{b}(z)} Z\_{i} Z\_{j} dz = 0 \quad , \qquad \qquad i \neq j \tag{49}$$

The orthonormal depth functions are obtained as

$$\overline{Z}\_{\rangle}(\mathbf{z}) = Z\_{\rangle} \wedge \sqrt{\left\{ \overline{\phi}^{(\prime)} \right\}^{T} \left[ \mathbf{K}^{\text{III}} \right] \left\{ \overline{\phi}^{(\prime)} \right\}} \tag{50}$$

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**

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

w

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

216 Environmental Applications of Remote Sensing

Eq. (47) satisfy the following orthogonality condition:

1 2

*D D*

1 1 0 , ( ) ( )

*i j i j D b*

*z z*

 r

<sup>1</sup> 0

The orthonormal depth functions are obtained as

r

eigenvectors {*ϕ*¯( *<sup>j</sup>*)

where

( ){ } { } 2 II I 2 III KK K *<sup>r</sup>*

( ) ( { }) <sup>2</sup> , , 1,2,.. *<sup>j</sup> rj k jn*

Reiss model that is based on Eqs. (30)–(32). Having obtained the eigenvalues *krj*

( ) { } { } ( ) ( ) *<sup>T</sup> <sup>j</sup> <sup>j</sup>* Ζ ψ *z z* = f

{ } ( ) 1 2 ( ), ( ), , ( ) *<sup>T</sup> <sup>n</sup>* ψ *z zz z* = yy

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 approx‐ imating 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 KIII -orthogonal, it can easily be shown that the eigenfunctions in

*Z Z dz Z Z dz i j*

{ } { } ( ) ( ) () / *<sup>T</sup> j j III Zz Z K j j* <sup>=</sup> f

+= ¹ ò ò (49)

 fé ù

It may be noted here that the eigenproblem in Eq. (45) remains linear unlike the Porter and

}, the eigenfunctions /depth modes may be written, using Eq. (40), as

 y

f

ff

é ùé ù é ù - = *<sup>k</sup>* ë ûë û ë û (45)

= (46)

(47)

L (48)

ë û (50)

2

and the

In our Rayleigh–Ritz model, the first task is to compute the symmetric matrices K<sup>I</sup> , KII , and KIII 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* = 2*n*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:


### **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.

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 redistrib‐ ution of power among the modes due to the presence of the rectangular hump.

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 Fig. 4. Transmission loss of a shallow-water column with a (a) rectangular hump on the seabed and (b) flat-bottom surface **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

among the modes due to the presence of the rectangular hump.

 Fig. 5. Spectra of modal efficiencies for (a) seabed with a rectangular bump and (b) a flat seabed **Figure 5.** Spectra of modal efficiencies for (a) seabed with a rectangular bump and (b) a flat seabed

Fig. 5. Spectra of modal efficiencies for (a) seabed with a rectangular bump and

#### **14.3. Down-sloping bottom** (a) (b)

geometry are shown in Fig. 6.

geometry are shown in Fig. 6.

**14.1. Isovelocity case**

218 Environmental Applications of Remote Sensing

**Figure 3.** Idealized ocean waveguide

especially on the right of the hump.

**14.2. Rectangular hump**

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

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

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 redistrib‐

ution of power among the modes due to the presence of the rectangular hump.

is chosen appropriately so that the modal error is less than 5%.

<H2>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 or down-slope, depending on the location of the source with respect to the slope. First we consider the down-sloping case where the 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 downsloping 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. (b) a flat seabed <H2>Down-Sloping Bottom

ocean-bottom slopes down from 100 m to 230 m over a distance of 600 m. The details of the

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 or 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

20

20

<H2>Down-Sloping Bottom

geometry are shown in Fig. 6.

among the modes due to the presence of the rectangular hump.

(b) flat-bottom surface

seabed

Fig. 4. Transmission loss of a shallow-water column with a (a) rectangular hump on the seabed and

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

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

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 or down-slope, depending on the location of the source with respect to the slope. First we consider the down-sloping case where the

 Fig. 6. Transmission loss of a shallow-water column with (a) down-sloping bottom and (b) flat bottom of depth 100 m **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. Fig. 7. Comparison of TL at 36 m depth of an ocean with down-sloping bottom and a flat bottom

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

Fig. 8. Modal power spectrum of shallow-water column with (a) down-sloping bottom and

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. 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.

#### <H2>Up-Sloping Bottom **14.4. Up-sloping bottom**

**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

20

Panel (a) shows the TL with the down-sloping bottom. Panel (b) shows the TL for a water column with the flat bottom at 100 m. Both results are for the source frequency of 150 Hz. Notice the

Fig. 6. Transmission loss of a shallow-water column with (a) down-sloping bottom and (b)

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

(a) (b)

Fig. 4. Transmission loss of a shallow-water column with a (a) rectangular hump on the seabed and

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

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

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 or 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

among the modes due to the presence of the rectangular hump.

(b) flat-bottom surface

<H2>Down-Sloping Bottom

geometry are shown in Fig. 6.

220 Environmental Applications of Remote Sensing

flat bottom of depth 100 m

seabed

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. 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.

21 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.

we studied. **Figure 9.** Transmission loss of shallow-water column with up-slope bottom at 105 Hz and 150 Hz

Fig. 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

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

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

#### <H2>Object in the Water Column **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. 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. (a) (b) Fig. 10. Modal spectrum of shallow-water column with up-sloping bottom (a) 105 Hz (b) 150 Hz <H2>Object in the Water Column

22

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

22

(a) (b)

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

<H2>Shallow-Water Column with Rippled Top Surface

conditions. In contrast, for the case of down slope the power flow is good for all the frequencies that

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

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.

conditions. In contrast, for the case of down slope the power flow is good for all the frequencies that

object, depending on the frequency of operation. This is because of the interference phenomena

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

we studied.

we studied.

<H2>Object in the Water Column

 Fig. 11. TL of a shallow-water column with an object (a) 135 Hz (b) 150 Hz **Figure 11.** TL of a shallow-water column with an object (a) 135 Hz (b) 150 Hz Fig. 11. TL of a shallow-water column with an object (a) 135 Hz (b) 150 Hz

#### <H2>Shallow-Water Column with Rippled Top Surface **14.6. Shallow-water column with rippled top surface** <H2>Shallow-Water Column with Rippled Top Surface

involving the object and the boundaries of the waveguide.

22

(a) (b)

22

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

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

conditions. In contrast, for the case of down slope the power flow is good for all the frequencies that

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

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

(a) (b)

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

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

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

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

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

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

interference phenomena involving the object and the boundaries of the waveguide.

we studied.

we studied.

<H2>Object in the Water Column

<H2>Object in the Water Column

**14.5. Object in the water column**

222 Environmental Applications of Remote Sensing

involving the object and the boundaries of the waveguide.

<H2>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. 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 wave‐ guide. 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. 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.

22 **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

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.

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.

Fig. 13. Sound-speed profile ( *S*1) used for our study

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

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

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

the source signal wavelength are key factors influencing the physics.

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

shown in Fig. 13.

224 Environmental Applications of Remote Sensing

**Figure 13.** Sound-speed profile (S1) used for our study

kind of sound-speed variation does not influence the TL much.

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

#### **14.8. Shallow-water column with depth-dependent sound speed and a rectangular bump on seabed** The TL for our geometry with the sound-speed profile given in Fig. 13 is shown in Fig. 14. The

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

Next, we consider the case of shallow-water ocean with depth-dependent sound speed and a rigid rectangular hump on the seabed. Next, we consider the case of shallow-water ocean with depth-dependent sound speed and a rigid 23 result for source frequency of 150 Hz is shown in Panel (a) and that corresponding to isovelocity is

rectangular hump on the seabed.

bump on the seabed, and (b) depth-dependent sound speed and flat bottom **Figure 15.** TL for the shallow-water ocean with (a) depth-dependent sound speed and a rectangular bump on the sea‐ bed, and (b) depth-dependent sound speed and flat bottom

Fig. 15. TL for the shallow-water ocean with (a) depth-dependent sound speed and a rectangular

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 depthdependent ocean. 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 depthdependent ocean.

Finally, we consider the case of shallow-water ocean with depth-dependent sound speed (Fig. 13)

<H2>Shallow-Water Column with Depth-Dependent Sound Speed and Rippled Top Surface

and a rippled air–water surface. The results are shown in Fig. 16.

(a) (b)

water interface, and (b) depth-dependent sound speed

24

Fig. 16. TL for the shallow-water ocean with (a) depth-dependent sound speed and a rippled air–

rectangular hump on the seabed.

#### **14.9. Shallow-water column with depth-dependent sound speed and rippled top surface** <H2>Shallow-Water Column with Depth-Dependent Sound Speed and Rippled Top Surface

Next, we consider the case of shallow-water ocean with depth-dependent sound speed and a rigid

Fig. 15. TL for the shallow-water ocean with (a) depth-dependent sound speed and a rectangular

We observe that, for the case when the source frequency is 60Hz, the presence of the small

bump on the seabed, and (b) depth-dependent sound speed and flat bottom

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. 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.

Fig. 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 **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. 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.

#### <H1>**Conclusion 15. Conclusion**

24 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 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.

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

<H2>Shallow-Water Column with Depth-Dependent Sound Speed and Rippled Top Surface

Next, we consider the case of shallow-water ocean with depth-dependent sound speed and a rigid

Fig. 15. TL for the shallow-water ocean with (a) depth-dependent sound speed and a rectangular

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-

Finally, we consider the case of shallow-water ocean with depth-dependent sound speed (Fig.

Finally, we consider the case of shallow-water ocean with depth-dependent sound speed (Fig. 13)

24

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

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.

Fig. 16. TL for the shallow-water ocean with (a) depth-dependent sound speed and a rippled air–

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

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

ripples can have a significant impact on the underwater propagation characteristics.

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

13) and a rippled air–water surface. The results are shown in Fig. 16.

(a) (b)

water interface, and (b) depth-dependent sound speed

and (b) depth-dependent sound speed

significant impact on the underwater propagation characteristics.

and a rippled air–water surface. The results are shown in Fig. 16.

bump on the seabed, and (b) depth-dependent sound speed and flat bottom

rectangular hump on the seabed.

226 Environmental Applications of Remote Sensing

dependent ocean.

<H1>**Conclusion**

**15. Conclusion**

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 finitedifference-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, *IR*(*pe*), from the radiation boundary of a finite element is represented by the second integral in that equation; i.e.,

$$I\_R(p\_\epsilon) = \frac{1}{2} \sum\_{m=1}^{M} \int\_{\mathcal{\rho}} \frac{1}{\rho} \alpha\_m p\_m^2 \, r d\boldsymbol{z} \tag{51}$$

where *M* denotes the number of propagating modes, *αm* the damper coefficient associated with the *m-*th mode [see Eq. (5)], *pm*(*z*) the pressure associated with the *m-*th normal mode, and *SRe* the element surface on the radiation boundary (see Fig. 1).

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

$$p\_m(\mathbf{z}) = a\_m f\_m(\mathbf{z}) \qquad m = 1, 2, \dots, M,\tag{52}$$

where *f <sup>m</sup>*(*z*) denotes the normal-mode function and *am* the modal coefficient. Using the finite element representation, the modal pressure on the radiation boundary may be written as

$$p\_m(\mathbf{z}) = \sum\_{\varepsilon} \begin{bmatrix} N \end{bmatrix}^T \{ \overline{p}\_{em} \}\_{\prime} \tag{53}$$

where *N* denotes the shape functions and {*p*¯*em*} 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\_c) = \frac{1}{2} \sum\_{m=1}^{M} \left\{ \overline{p}\_{em} \right\}^T \left\{ \overline{R}\_{cm} \right\} \left\{ p\_{cm} \right\}, \tag{54}$$

where

$$
\left[\overline{R}\_{eu}\right] = \alpha\_{\
u} \int\_{\mathcal{S}\_{\mu}} \frac{1}{\rho} \left[\boldsymbol{N}\right]^{\intercal} \left[\boldsymbol{N}\right] d\boldsymbol{S}.\tag{55}
$$

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

$$\sigma\left\{\overline{p}\_{um}\right\} = a\_m \left\{f\_{zm}(\mathbf{z}\_1), f\_{zm}(\mathbf{z}\_2), \dots, f\_{zm}(\mathbf{z}\_n)\right\}^T = a\_m \left\{f\_{zm}\right\},\tag{56}$$

where *f zm*(*zj* ) 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} \sum\_{m=1}^{M} a\_m^2 R\_{cm'} \tag{57}$$

where

$$\mathcal{R}\_{cm} = \left\{ f\_{zm} \right\}^T \left[ \overline{\mathcal{R}}\_{cm} \right] \left\{ f\_{zm} \right\}. \tag{58}$$

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

$$
\mu\_{\rangle} = a\_{\rangle} \sin(k\_{z\rangle} z) \tag{59}
$$

where

where *N* denotes the shape functions and {*p*¯*em*} 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

{ } { }{ } <sup>1</sup>

<sup>1</sup> .

*T*

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

with the radiation boundary. Now, using Eqs. (55) and (56), the functional in Eq. (54) may be

<sup>=</sup> å (54)

é ùé ù ë û ò ë ûë û (55)

<sup>=</sup> å (57)

ë û (58)

*em m zm zm zm n m zm p afzfz fz af* = = K (56)

) denotes the *j-*th nodal value of the *m-*th eigenmode on a finite element in contact

,

2 1 ( ) , *<sup>M</sup> <sup>T</sup> R e em em em m Ip p R p* =

*em m*

 é ù = a

*Re*

{ } ( ), ( ), , ( ) 1 2 { }, *<sup>T</sup>*

1 2 2 1

{ } { }. *<sup>T</sup> R f Rf em zm em zm* = é ù

**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

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

*M R m em m I aR* =

*S R N N dS* r

see Eq. (15)]

228 Environmental Applications of Remote Sensing

where

written as

where *f zm*(*zj*

written as

where

$$k\_{zj} = \frac{\left(j - 0.5\right)\pi}{D}, j = 1, 2, \dots \tag{60}$$

where *aj* is chosen to normalize the mode functions.

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

$$\begin{aligned} \, \nu\_{/}(z) &= a\_{/} \sin k\_{z/} z & 0 \le z \le D\_{1} \\ &= a\_{/} \frac{\sin k\_{z/} D\_{1}}{\cos k\_{b \nearrow} (D\_{2} - D\_{1})} \cos k\_{b \nearrow} (D\_{2} - z) & D\_{1} \le z \le D\_{2} \end{aligned} \tag{61}$$

where *kzj* <sup>2</sup> <sup>=</sup>*<sup>ω</sup>* <sup>2</sup> / *<sup>c</sup>* <sup>2</sup> <sup>−</sup>*krj* 2 , *kbzj* <sup>2</sup> <sup>=</sup>*<sup>ω</sup>* <sup>2</sup> / *cb* <sup>2</sup> <sup>−</sup>*krj* 2 and *krj* is the solution of the transcendental equation

$$\frac{k\_{z\bar{\jmath}}\rho\_b}{k\_{bz\bar{\jmath}}\rho} = \tan k\_{z\bar{\jmath}} D\_1 \tan k\_{bz\bar{\jmath}} (D\_2 - D\_1) \tag{62}$$

**Figure 17.** A bounded two-layer waveguide

In the case of a Pekeris waveguide (for which *D*2→*∞* in Fig. 17), the trial functions are given by

$$\begin{aligned} \, \vert \boldsymbol{\nu}\_{/}(\boldsymbol{z}) = \boldsymbol{a}\_{/} \sin k\_{\boldsymbol{z}/} \boldsymbol{z} & \quad \quad 0 \le \boldsymbol{z} \le \boldsymbol{D}\_{1} \\ &= \boldsymbol{a}\_{/} \sin \langle k\_{\boldsymbol{b}\boldsymbol{z}} \boldsymbol{D}\_{1} \rangle e^{-\stackrel{\scriptstyle \cdot \vert k\_{\boldsymbol{b}\boldsymbol{z}} \cdot \left(\boldsymbol{D}\_{1} - \boldsymbol{z}\right)}} \quad \boldsymbol{D}\_{1} \le \boldsymbol{z} \le \boldsymbol{\infty} \end{aligned} \tag{63}$$

where *kzj* <sup>=</sup> *<sup>ω</sup>* <sup>2</sup> / *<sup>c</sup>* <sup>2</sup> <sup>−</sup>*krj* 2 , *kbzj* =*i krj* <sup>2</sup> <sup>−</sup>*<sup>ω</sup>* <sup>2</sup> / *cb* 2 and *krj* is the solution of the transcendental equation

$$\tan k\_{z\rangle} D\_1 = -\frac{\mathrm{i}\,\rho\_b k\_{z\rangle}}{\rho k\_{bz\rangle}} \tag{64}$$

Note that discrete guided modes in the water column exist only for the case when *krj* <sup>2</sup> <sup>&</sup>gt;*<sup>ω</sup>* <sup>2</sup> / *cb* 2 .

There are two cases to consider:

Case 1: *k* <*kb*

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* >*kb*

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 men‐ tioned 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* <*kb*

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

Case 2: *k* >*kb*

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.

#### **Acknowledgements**

In the case of a Pekeris waveguide (for which *D*2→*∞* in Fig. 17), the trial functions are given

1

*ik D z*


*a kDe D z*

( ) sin 0 sin( ) *bzj*

*z a kz z D*

<sup>1</sup> tan *b zj zj*

Note that discrete guided modes in the water column exist only for the case when *krj*

*k D*

*j j zj*

y

230 Environmental Applications of Remote Sensing

2

, *kbzj* =*i krj*

sound signal over long distances in the water column.

attention to discrete guided modes as shown above.

boundary. Therefore, the underlying physical processes are different.

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

long range propagation, the modes of significance are the discrete guided modes.

*j bzj*

<sup>2</sup> <sup>−</sup>*<sup>ω</sup>* <sup>2</sup> / *cb* 2

( ) 1 1

*bzj*

*i k*

*k* r

r

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.

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

Since our interest is in long-range sound transmission in the water column, we have restricted

One should observe that our two-layer waveguide problem does not share the above men‐ tioned behavior. Note that the two-layer waveguide is terminated at the bottom by a rigid

In this case, the spectrum consists of discrete guided modes and surface modes. However, for

1

<sup>=</sup> £ £¥ (63)

and *krj* is the solution of the transcendental equation

= - (64)

<sup>2</sup> <sup>&</sup>gt;*<sup>ω</sup>* <sup>2</sup> / *cb* 2 .

by

where *kzj* <sup>=</sup> *<sup>ω</sup>* <sup>2</sup> / *<sup>c</sup>* <sup>2</sup> <sup>−</sup>*krj*

There are two cases to consider:

Case 1: *k* <*kb*

Case 2: *k* >*kb*

Case 1: *k* <*kb*

Case 2: *k* >*kb*

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.

#### **Author details**

Saba Mudaliar1\*, C.P. Vendhan2 and C. Prabavathi3

\*Address all correspondence to: saba.mudaliar@us.af.mil

1 Sensors Directorate, Air Force Research Laboratory, WPAFB Ohio, USA

2 Indian Institute of Technology Madras, Chennai, India

3 Independent researcher, Dayton, Ohio, USA

#### **References**


[25] Antoine, X., H. Barucq, and A. Bendali (1999) Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape, *J. Math. Anal. Appl*., Vol. 229, pp. 184–211.

[10] Evans, R.B. (1983) A coupled mode solution for acoustic propagation in a waveguide with stepwise depth variations of a penetrable bottom, *J. Acous. Soc. Am.*, Vol. 74, pp. 188–

[11] Tappert, F.D (1977) The parabolic approximation method, wave propagation and underwater acoustics, in *Wave Propagation and Underwater Acoustics*, J.B. Keller and J.S.

[12] DiNapoli, F.R. and R.L. Deavenport (1980) Theoretical and numerical Green's function field solution in a plane multilayered medium, *J. Acous. Soc. Am.*, Vol. 67, pp. 92–105.

[13] Zienkiewicz, O.C. and R.L. Taylor (1989) *The Finite Element Method*, *Vol. 1- Basic Formula‐*

[14] Bathe, K.J. (1982) *Finite Element Procedures in Engineering Analysis*, Prentice Hall, Inc.,

[15] Cook R.D., D.S. Malkus, M.E. Plesha, and R.J. Witt (2002) *Concepts and Applications of Finite*

[16] Schmidt, H. and F.B. Jensen (1985) A full wave solution for propagation in multilay‐ ered viscoelastic media with application to Gaussian beam reflection and fluid-solid

[17] Thompson, D.J. and N.R. Chapman (1983) A wide-angle split-step algorithm for the

[18] Lee, D., G. Botseas, and J.S. Papadakis (1981) Finite difference solution to the parabolic

[19] Huang, D. (1988) Finite element solution to the parabolic wave equation, *J. Acous. Soc.*

[20] Evans, R.B. and K.E. Gilbert (1985) Acoustic propagation in a reflecting ocean wave‐

[21] Dougalis, V.A., N.A. Kampanis, and M.I. Taroudakis (1998) Comparison of finite element and coupled mode solutions of the Helmholtz equation in underwater acoustics, *Proc. 4th European Conference in Underwater Acoustics*, eds. A. Alippi and G.B. Cannelli, CNR-

[22] Mitsoidis, D.A., N.A. Kampanis, and V.A. Dougalis (2008) Coupled mode and finite element approximations of underwater sound propagation problems in general

[23] Murphy, J.E. and S.A. Chin-Bing (1988) A finite element model for ocean acoustic

[24] Murphy, J.E. and S.A. Chin-Bing (1989) A finite-element model for ocean acoustic

propagation and scattering, *J. Acous. Soc. Am.*, Vol. 86, pp. 1478–1481.

guide with an irregular interface, *Comp. Math. Appls*. Vol. 11, pp. 795–805.

Papadakis, eds., Springer-Verlag, Berlin, pp. 224–287.

*tion and Linear Problems*, McGraw-Hill, New York.

*Element Analysis*, 4th ed., John Wiley, New York.

interfaces, *J. Acous. Soc. Am.*, Vol. 77, pp. 813–25.

parabolic equation, *J. Acous. Soc. Am.*, Vol. 74, pp. 1848–1854.

stratified environments, *J. Comput. Acous.*, Vol. 16, pp. 83–116.

propagation, *Math. Comput. Modelling*, Vol. 11, pp. 70–74.

wave equation, *J. Acous. Soc. Am.*, Vol. 70, pp. 795–800.

Englewood Cliffs, New Jersey.

*Am.*, Vol. 84, pp. 1405–1413.

IDAC, Rome, Vol. II, pp. 649–654.

195.

232 Environmental Applications of Remote Sensing

