Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research Open Access

Use of analytical and numerical scattering models to predict acoustic backscatter from aquatic organisms continues to increase in both fisheries management and ecosystem research. These models have been used to predict backscatter for classification of volume backscatter and target strength (TS; dB re 1 m²) of individual scatterers, and to corroborate acoustic measurements with directed net samples. The effectiveness of a model is generally evaluated for a specific species or survey, which makes selecting a model or set of models difficult for scientists who require accurate number density and abundance estimates for a suite of aquatic organisms. In this paper, models are defined as the analytical and numerical mathematical expressions implemented using computer algorithms to predict acoustic backscatter; predictions are the resulting backscatter amplitudes; and the digital representations of the target shape and properties (e.g., organism anatomy and morphometry, material properties, and boundary conditions) are used as input to the models.

Early observations that echoes from aquatic organisms can dominate other reflectors of underwater sound motivated the need to quantify and predict acoustic backscatter from biological targets (Midttun, 1984). Backscatter can be measured from biological targets in their natural environment (i.e., in situ) or under controlled experimental conditions (i.e., ex situ). Relationships of these measurements to animal lengths derived from trawl catches may be estimated empirically using statistical regressions, theoretically by comparisons with analytical and numerical models, or both (Medwin and Clay, 1998). Empirical approaches can be used to investigate the magnitude and relative importance of biological and physical factors that influence backscatter intensities, but a complete understanding of sound scattering by aquatic organisms requires a combination of empirical and theoretical methods (Henderson and Horne, 2007).

The use of models to predict backscatter by aquatic organisms rapidly expanded in the 1960s (Haslett, 1965) and numerical techniques to solve analytical models followed as computing capabilities improved (Francis and Foote, 2003). Applications of backscatter predictions to estimate fish and zooplankton abundance (Holliday, 1972; Lavery et al., 2007) have paralleled the development of anatomically and morphologically accurate representations of targets. Such representation of aquatic organisms has evolved from simple geometric shapes such as spheres, cylinders, and prolate spheroids (Anderson, 1950; Jech et al., 1995), to more accurate representations of fish anatomy (Clay and Horne, 1994; Reeder et al., 2004), zooplankton anatomy (Chu et al., 1993), and cephalopod anatomy (Lee et al., 2012). The availability of multifrequency acoustic data has facilitated an increase in use of theoretical models to separate echoes from fish with gas or oil inclusions (e.g., swimbladders) from fish without swimbladders, and zooplankton (Kloser and Horne, 2003). Multifrequency data continue to be used to separate echoes from multiple species within a survey region (Anderson et al., 2007; De Robertis et al., 2010; Woillez et al., 2012).

Choosing from among theoretical models those that are optimal for predicting backscatter from fish and zooplankton is not trivial. There are a number of models available, each with advantages and constraints. In this study, models were chosen based on historical usage and potential application to biological targets with high aspect ratios. The predictions of backscatter from these models have not been systematically compared using a common set of target shapes and properties.

In this paper, ten analytical and numerical models were applied to four target shapes with homogeneous composition and their predictions compared among models and to benchmarks (i.e., predictions of acoustic backscatter from exact or approximate analytical models) (Table I). These shapes and compositions were selected for comparison because they have exact analytical solutions or the approximations are known and accepted. In this paper, all models were implemented in computer algorithms and were numerically evaluated to within a predetermined precision (Sec. IV C). The shapes were a sphere, a sphere with a fluid shell, a prolate spheroid, and a finite cylinder. For each target shape, four target boundary conditions were modeled: fixed-rigid, pressure-release, gas-filled (air-filled), and weakly scattering. Material properties of the gas-filled target were selected to represent a gas-filled swimbladder and the properties of the weakly scattering target were selected to represent macrozooplankton (e.g., krill). The objectives of this study were to (1) predict TS of the simple targets as a function of frequency and insonification angle, (2) compare predictions from different models using the same target, and (3) identify advantages and constraints of each model.

Where possible, solutions of exact analytical models were used as benchmarks (Table I). This includes the geometries of the sphere and prolate spheroid where exact modal series solutions (MSSs) to the wave equation were used. In cases where an exact solution did not exist, either an approximate analytical solution was used as the benchmark or no benchmark predictions were made. The modal series–based deformed cylinder solution was used in this former case, as its limitations have been documented in earlier studies.

The exact modal series model used to describe the scattering by a sphere is obtained using separation of variables. The general solution for the scattered wave by a sphere can be expressed as

where r is the range (m), k_w is the acoustic wavenumber (k = 2π/λ) in the surrounding water, subscript w denotes the surrounding water, λ is wavelength (m), p_o is the incident pressure, and the value of the coefficient $An$ depends on the boundary conditions (Anderson, 1950). Note that our definition of A_n (see below) and Eq. (1) differs from Anderson (1950). P_n is the Legendre polynomial of degree n, and $hn1$ is the spherical Hankel function of the first kind, or the spherical Bessel function of the third kind, of order n. For backscattering, the scattering angle is $θ=π$ and the Legendre polynomial becomes P_n(−1) = (−1)ⁿ. For the far field, k_wr ≫ 1, $hn1(kwr)kwr≫1→[(−i)n+1/(kwr)]eikwr$ and the backscattered acoustic pressure can be expressed in terms of the scattering amplitude

where

is the backscattering form function. For these calculations, k_wa + 20 (rounded to the nearest integer) was used as the maximum number of terms for n, which provided the sufficient number of terms for convergence relative to the precision of 0.1 dB required for these comparisons (see Sec. IV C). Note that the scattering amplitude has dimension of length and is also known as scattering length (Medwin and Clay, 1998). TS is related to $fbs$ through the expression:

Note that the coefficient A_n given in Eq. (3) differs from that given in Anderson (1950) by a factor of −p_o(−i)ⁿ(2n + 1) due to the difference in our definition of A_n. However, our results are identical to Anderson.

For fixed-rigid spheres, the boundary condition requires that the normal velocity is zero at the water/sphere interface, which leads to

where j_n is the spherical Bessel function of the first kind, the prime denotes the derivative with respect to the argument, and a is the radius of the sphere.

For pressure-release spheres, the boundary condition requires zero acoustic pressure at the water/sphere interface, which leads to

(Rayleigh, 1945; Morse and Ingard, 1968).

For fluid-filled spheres, the acoustic pressure and the normal velocity are both non-zero at the water/sphere interface. Fluids are liquids and gases that do not support a shear wave. The coefficient $An$ (Anderson, 1950) is

where

$k1$ is the acoustic wavenumber inside the fluid-filled sphere, g is the density (ρ, kg m³) contrast $(ρt/ρw)$ (subscript t denotes the target), h is the sound speed (c, m s⁻¹) contrast $(ct/cw)$ of the sphere to the surrounding water, and y is the spherical Neumann function. In addition, letting gh→∞ and gh→0 for the fluid-filled sphere, Eq. (8) reduce to Eqs. (5) and (6) corresponding to the fixed-rigid and pressure-release cases, respectively.

For a gas-filled sphere (i.e., air-filled bubble), the density contrast of the fluid-filled sphere is near zero (e.g., g ≈ 0.0012 at the surface) and the sound speed contrast is less than unity (e.g., h ≈ 0.22 at the surface). Equations (7) and (8) were used with these g and h parameter values to compute the exact scattering amplitude.

For the weakly scattering sphere, the density and sound speed contrasts are both near unity. Equations (7) and (8) were used to compute the exact scattering amplitude.

For the spherical fluid shell with fluid interior, the acoustic pressure and normal velocity are finite at both the water/outer-sphere and the inner-sphere interfaces (Jones et al., 2009). In the benchmark case here, the 4 × 4 matrices in Jones et al. (2009) were transformed to 3 × 3 matrices to simplify benchmark computations by removing the unused coefficients for the inner-layer pressure (i.e., eliminating the 3rd columns of the original 4 × 4 matrices). The coefficient A_n can be expressed as

where all elements in the two 3 × 3 matrices are listed in the  Appendix.

The fixed-rigid spherical shell is acoustically equivalent to the fixed-rigid solid sphere and therefore a benchmark was not calculated for this case.

For the spherical fluid shell with pressure-release interior, the density contrast between the interior fluid and shell is zero (i.e., g₃₂ = 0). A_n can be expressed as

where elements in the two 2 × 2 matrices are listed in the  Appendix.

For the spherical fluid shell with gas interior, the density and sound speed contrasts of the gas relative to the shell are both much smaller than unity. Equations (3) and (9) were used to calculate TS.

In this case, both the fluid shell and the inner fluid sphere are weakly scattering. The resulting density and sound speed contrasts—$g21$⁠, $g32$⁠, $h21$⁠, and $h32$—are all near unity (1.002). Equations (3) and (9) were used to calculate TS.

The prolate spheroid modal series (PSMS) model was used as the benchmark for the fixed-rigid, pressure-release, and gas-filled prolate spheroid targets.

Solving the scalar wave equation in spheroidal coordinates (ξ, η, ϕ) and applying the boundary condition determines the scattering amplitude of the scattered wave from a prolate spheroid with a major radius a and a minor radius b (Flammer, 1957; Silbiger, 1963; Yeh, 1967; Skudrzyk, 1971). A spheroid surface is given by ξ = ξ_w = constant and the relationships among a, b, and ξ are a = ξ_wq and $ξw=1/1−(b/a)2$⁠, where q is the semi-focal-length of the spheroid.

The scattering amplitude (f_sc) is expressed by the prolate-spheroidal wave function approximated in the far field as

where $hw≅kwq$⁠, 2q is the distance between focal points, (⁠$θ,ϕ$⁠) and (⁠$θ′,ϕ′$⁠) are the spherical angle coordinates of the scattered and incident wave, $εm$ is the Neumann factor, S_mn is the angular spheroidal wave function of the first kind of order m and degree n, and N_mn is its normalization factor. For backscatter, $θ=π−θ′$ and $ϕ=π+ϕ′$⁠. Values of the unknown coefficients A_mn are determined by the boundary conditions required for both the pressure field and the normal component of displacement across the scattering surface. The maximum values of m and n in Eq. (11) were estimated by Furusawa (1988) as

In the backscattering case, θ′ = π−θ, ϕ′ = π, and ϕ = 0 in Eq. (11). The TS of a scatterer is defined in Eq. (4) using f_bs = f_sc(θ,0|π−θ,π).

Only a small number of parameters are required for the calculations. For example, four parameters (k_w, q, ξ_w, θ) are required for the pressure-release and rigid fixed model calculations for the backscattering case and another two parameters (k₁,g) are required for the fluid case.

The PSMS is correct for all values of h_ps and for all angles of incidence naturally encountered. However, computer precision limits the maximum value of kq to approximately 12. In the gas-spheroid model, the gas parameters were applied to the fluid model and calculated over a frequency range of f = 12 to 40 kHz. The pressure-release model was calculated over a frequency range of f = 12 to 80 kHz. Comparisons between the model and experimental results have shown good agreement (Sawada et al., 1997). Comparisons between the PSMS and the deformed cylinder model (DCM) with vacant, rigid, and fluid conditions and between the PSMS and the boundary element method (BEM) agreed well (Ye et al., 1997; Okumura et al., 2003).

Spence and Granger (1951) and Senior (1960) showed that coefficients for the rigid spheroids can be expressed as

where $Rmn(i)$ is the radial spheroidal wave function of the ith kind.

Senior (1960) showed that the coefficients for the pressure-release spheroid can be expressed as

Coefficients for the fluid-filled spheroid were solved using the simultaneous equation

where m = 0,1,2,…, Σ′ indicates a summation with respect to the same parity of n and l, and $Knlm(i)$ is shown as

and

where $Rmn′=∂Rmn/∂ξ$⁠. In the case that $ht≅hw$⁠, A_mn can be further simplified to

which was used for calculations of the benchmark.

The PSMS model did not converge for a gas spheroid at “high” ka values, thus no benchmark model was used for this target.

For the weakly scattering spheroid, the g and h values are near unity. Equations for the fluid-filled prolate spheroid described in Sec. II A 2 c were used to compute scattering amplitudes.

There is no exact analytical model for a straight finite cylinder. For this target, the modal series–based DCM was used as the benchmark for the homogeneous rigid, pressure-release, gas-filled, and weakly scattering cylinders at broadside incidence and as the benchmark for these targets at f = 38 kHz from θ = 70° (20° off broadside) to 90° (broadside incidence).

The modal series–based DCM was developed in Stanton (1988, 1989). The solution involves integrating the scattered pressure per unit length along the length of the deformed cylinder. For the case of a straight finite cylinder, the integral reduces to the expression for the backscattering amplitude, f_bs:

where L is the cylinder length, k is the wavenumber in the surrounding water, and θ is the angle between the straight cylinder axis and the incident plane wave [Eq. (28) in Stanton (1988) or Eq. (12) in Stanton (1989)]. Here the coefficient B_m are determined by matching interior and exterior solutions to the wave equation at the boundary of the scatterer. These coefficients are specific to the type of material property of the scatterer (pressure-release, fixed-rigid, and fluid/gas) relative to the surrounding fluid and are given below for each of those cases.

The modal series–based DCM solution assumes that (1) end effects are negligible and (2) the direction of the tangent to the cylinder axis, the cross-sectional radius, and the material properties change slowly with respect to position along the cylinder axis. The first assumption restricts this use of the DCM to geometries where the direction of the incident sound wave and scattering is normal or near normal to the tangent of the cylinder axis, and the cylinder is elongated where the aspect ratio is much greater than unity. For applications to the straight undeformed cylinders of homogeneous properties, the second assumption is automatically met. Note that the DCM was later generalized for any kernel (such as modal series–based, ray-based, DWBA-based) in Stanton (1992) and, for the DWBA-based-DCM applications described below, the solution is not restricted to the above assumptions.

The backscattering amplitude of a rigid and fixed straight cylinder is [Eq. (34), Stanton, 1988]

where ε_o = 1 and ε_m = 2 for m = 1, 2, 3, 4…, J_m(x) is the cylindrical Bessel functions of order m of the first kind, the prime represents the derivative with respect to the arguments, and $K=k sin θ$⁠, H_m is the cylindrical Bessel function of the third kind (Hankel function). TS is computed using Eqs. (4) and (19).

The backscattering amplitude of a pressure-release finite cylinder is [text after Eq. (17) of Stanton, (1988)]:

TS values were computed using Eqs. (4) and (20).

From Stanton (1988), the coefficients B_m in Eq. (18) of this paper for this fluid-fluid boundary condition are

where

N_m is the cylindrical Bessel function of the second kind of order m, and $K′=K/h$⁠, TS values were computed using Eqs. (4), (18), (21), and (22) and the density and sound speed values of gas from Sec. IV B. Note that acoustically, gas is a fluid as it does not support a shear wave.

The fluid-fluid boundary conditions here are identical to those of the gas-filled cylinder. TS values were computed using Eqs. (4), (18), (21), and (22) and the material properties of the weakly scattering target given in Sec. IV B.

The acoustic BEM uses an integral form of the Helmholtz wave equation, in which the acoustic pressure p at any point is expressed in terms of the pressure and normal displacement fields on the surface S of the scatterer (Chen and Scheikert, 1963; Chertock, 1964; Copley, 1967). By evaluating this integral equation at each node of a mesh of elements that collectively span S, a system of simultaneous equations for the pressures and displacements at those nodes is produced. For the pressure-release and rigid surfaces, the pressure and normal displacement, respectively, are known, and the system of equations may be solved for the other variable. Where S encloses a second fluid, the interior form of the standard integral equation provides a second system of equations, which can be solved simultaneously with the first for both surface pressure and normal displacement. This solution may then be used in the original integral equation to determine the pressure at any other point.

The particular implementation in this study used the approach of Burton and Miller (1971) to overcome the known problem of singularities in the exterior equation at certain critical frequencies (Copley, 1968; Pierce, 1992), in which the standard integral equation has a multiple of its normal derivative at the surface added to the integral. Further details of the method can be found in Francis (1993) and Francis and Foote (2003).

In this formulation, integrals are evaluated numerically by Gaussian quadrature over six-noded triangular and eight-noded quadrilateral elements that together span S; on each element the acoustic and geometric variables are expressed in terms of nodal values using quadratic interpolation (Zienkiewicz and Taylor, 1989). For this interpolation to closely represent the acoustic waveforms, it is recommended that adjacent nodes should be no more than $16$ of a wavelength apart (Bonnet, 1995).

This condition requires finer meshes as the frequency increases. For example, for an ellipsoid at 400 kHz, a triangular mesh covering the whole surface would need to consist of over 18 000 elements and 37 000 distinct nodes. To solve the dense, complex system of equations for the surface pressures at this number of nodes would present a formidable computing challenge. The size of the problem may be reduced by exploiting mirror symmetry in the plane containing the ellipsoid axis and the incident wave direction. To achieve further reduction, a superposition of solutions, with symmetry and anti-symmetry, has been implemented. For each of the other two principal planes, the problem is rendered symmetric by adding the mirror image of the incident wave, and then anti-symmetric by reversing the sign of this mirrored wave. The four cases that arise are solved separately and the resulting solutions are superimposed to solve the original problem using an incident wave of four times the amplitude. As a result, only one octant of the ellipsoid (or cylinder) needs to be meshed, requiring less than 5000 nodes at 400 kHz. This procedure was verified by comparing its results at a low frequency (12 kHz) with those obtained using a corresponding mesh of the complete surface, i.e., one derived by reflecting the octant mesh in the three coordinate planes in turn. Agreement (to the eighth significant figure) in the values of the scattered pressure was confirmed.

Tests for convergence were conducted at f = 12, 18, 38, 70, 120, 200, and 400 kHz for the rigid, pressure-release, and weakly scattering cases of the prolate spheroid and cylinder, using series of meshes with decreasing nodal separation. For the rigid and pressure-release targets, convergence was achieved to within margins of ± 0.02 dB at all frequencies, except at 200 kHz for the rigid cylinder, where the margins were ± 0.2 dB. Convergence was less consistent for the weakly scattering targets. Margins of ±0.1 dB were achieved for the spheroid at all frequencies up to and including 70 kHz; thereafter the margins were ±0.3, ±1.5 and ±3.5 dB at 120, 200, and 400 kHz, respectively. For the cylinder, margins of ±0.1 dB were achieved at all frequencies except 120 kHz (±2 dB); it may be noted that this is close to a null in the TS response, where the subtraction of two nearly equal matrices in the BEM makes the computations more sensitive to error.

The FEM was used to solve the inhomogeneous Helmholtz equation in the frequency domain to give the steady-state acoustic pressure throughout the simulation volume (Ihlenburg, 1998). The scattering objects were surrounded by a homogenous layer, and then by a perfectly matched layer (PML), which acted as an absorbing boundary (Berenger, 1996; Zampolli et al., 2007). The far-field backscattered pressure was estimated by evaluating the Helmholtz-Kirchhoff integral (Pierce, 1989) on the boundary between the homogenous layer and PML. This method is able to simulate the scattering from arbitrarily complex fluid-like inhomogeneous objects.

The fixed-rigid and pressure-release targets were modeled as a void in the model domain with appropriate boundary conditions. TS predictions for the gas-filled and weakly scattering targets were obtained by varying the model material properties to form the scattering object. The homogenous layer surrounding the scattering object was at least one-half wavelength thick, while the PML was one wavelength thick. The FEM elements were Lagrangian of quadratic order. The mesh density was set to give at least ten nodes per wavelength so as to adequately represent the acoustic wave (e.g., Ihlenburg, 1998). Tests for convergence were carried out at 38 kHz for the rigid, pressure-release, and weakly scattering cases of the sphere. Convergence of the solution occurred with ten or fewer nodes per wavelength for all three sphere cases. The ten-node per wavelength mesh density was maintained for all models. Symmetry in the model and scattered pressure field was used to reduce the computation load where possible.

The FEM implementation was provided by the Acoustics Module of the commercially available COMSOL Multiphysics software version 3.5 (COMSOL, 2008). The methods were implemented using the “time-harmonic scattered wave” option where the incident wave is used to generate an equivalent radiation condition on the boundary of the fixed-rigid and pressure-release scattering targets. The gas-filled and weakly scattering targets were simulated using a two domain approach where the total pressure field inside the target was calculated and coupled to a scattered field in the surrounding homogenous layer and PML. The coupling was achieved by setting the pressure on the boundary of the object equal to the total pressure in the homogenous layer, and ensuring the continuity of the normal derivative of velocity between the two domains.

The FMM uses a two-dimensional (2-D) conformal mapping approach to describe scattering by axisymmetric, irregular, finite-length bodies of revolution (Reeder and Stanton, 2004). The model conformally maps the coordinate variables of the original coordinate system to a new orthogonal coordinate system using the Newton–Raphson method, where the new radial coordinate exactly coincides with the target surface. The solutions to the transformed Helmholtz equation are a general solution for the total pressure in the case of far-field scattering by a finite body of revolution. This model has been shown to be accurate (i.e., within 1 dB of the exact solution) when predicting scattering by smooth, symmetric bodies over a wide range of frequencies (i.e., Rayleigh through geometric scattering regions), angles (monostatic and bistatic), aspect ratios, and boundary conditions (Reeder and Stanton, 2004).

The conformal mapping function has the form

where c_n are the conformal mapping coefficients, and ρ is related to the radial and polar angular variables in the new coordinate system. The scattering amplitude in the new coordinate system is

where b_nm are the far-field scattering coefficients that depend on the shape and material properties of the scatterer, $Pnm$ is the associated Legendre function, r(u,w) is the new radial coordinate, and g(u,w) is a function of the new system. The far-field scattering coefficients, b_nm, are determined after solving the transformed Helmholtz equation and satisfying the boundary conditions at the surface of the target. The values of m and n differed among targets, but in all cases exceeded ka + 3. Details and definitions of the parameters are in Reeder and Stanton (2004).

Given the eccentricity of shape and high frequencies at which these TS values are computed, approximations were required at the higher frequencies to obtain a stable result.

The KA is a high-frequency model in which diffraction is neglected. The pressure field on the scattering surface S is assumed to be known a priori. For a pressure-release surface, where the acoustic pressure is zero, the normal component of particle velocity on S is assumed to be equal to that of the incident field on the directly insonified part of S, and zero on the geometrically shadowed part. Similarly, for a rigid surface, where the normal velocity is zero, the pressure on S is assumed to be that of the incident wave on the insonified part and zero on the shadowed part of S. The scattered pressure at any point exterior to S can then be obtained from the Helmholtz integral equation. The form that this integral takes for the far-field backscattering amplitude may be found in Foote (1985) and Foote and Francis (2002).

In the model used here, the integration is performed numerically using the same procedure as in the BEM: Gaussian quadrature over curvilinear surface elements, where the position vector and the associated normal at each Gaussian point are calculated by quadratic interpolation from the nodal coordinates. Accuracy of the numerical integration depends on the order of the polynomial fit assumed in the Gaussian quadrature and on the distance between adjacent nodes. With Gaussian quadrature of order 2 or higher, a good representation of the waveform should be obtained for a nodal separation up to one-sixth of a wavelength. In the most extreme case considered here – the spheroid at 400 kHz – this condition leads to meshes of over 18 000 triangular elements. The KA does not need to solve for the pressure field on the surface, and requires an evaluation of a single integral, which significantly reduces computational requirements relative to the BEM.

The KRM calculates backscatter as a function of target length, acoustic frequency, and angle of incidence (Clay and Horne, 1994). Backscatter is calculated from targets using a low-mode cylinder solution for frequencies near resonance of the swimbladder and a Kirchhoff-ray approximation at high frequencies from digital representations of the target bodies and inclusions. The Kirchhoff-ray approximation accounts for local reflection coefficients of front and back interfaces of penetrable scatterers. Inclusions include gas- or lipid-filled swimbladders in fish, oil droplets in zooplankton, and skeletal elements. The digital target morphology was used to construct cylindrical volume elements. Digital resolution and element spacing in the KRM model were fixed at 1 mm for all targets and frequencies, as that is what is done operationally. For the higher frequencies, such as 400 kHz, this spacing is larger than $16$ of a wavelength which may lead to errors. Backscatter from each morphological feature, both in the low and high frequency regions is computed as a finite cylinder using the deformed finite cylinder model of Stanton (1989).

For the gas-filled targets with ka values ≤0.2, where a is the equivalent radius of the cylindrical elements, a low-mode (M = 0) cylinder solution was used (Clay, 1991, 1992). For the gas-filled targets and ka > 0.2, a Kirchhoff-ray approximation was used, which is modified from Clay (1992). For the weakly scattering targets at all ka, a Kirchhoff-ray approximation was used.

For angles from 0° (end-on incidence) to 69°, non-benchmark TS was calculated for the cylinder using the modal series–based DCM. The DCM was also used for TS predictions of the prolate spheroid targets. Since this involves cylinders deformed beyond the shape of the straight cylinders described in Eq. (18), integration of the scattering per unit length was performed along the length of the target, using Eq. (8) of Stanton (1989).

In the method used here, the integration is performed numerically using digitized shapes of the benchmark targets. Digitizing includes (i) digitally slicing elongated targets cross-wise and at closely spaced points along the major axis to form thin disk-like volume elements and (ii) defining their radius, position along the axis and the orientation of their axis. The digital resolution and the element length along the axis were approximately 0.3 mm for cylinders and 0.6 mm for the prolate spheroids. Here, the element lengths are less than λ/6 for all frequencies.

The DWBA, an approach that had long been used in other fields such as quantum mechanics, was first applied to the case of weakly scattering zooplankton by Chu, Stanton, and colleagues in a series of papers beginning in 1993 (Chu et al., 1993; Stanton et al., 1993). In this particular application, the “distortion” involved accounting for the change in wavenumber of the incident field within the body due to the material properties of the body. The DWBA is applicable only to weakly scattering targets where the material properties of the target are to within approximately 5% of the surrounding fluid and where there are no shear waves. The distinct advantage of this approach is that it is not restricted to scattering geometry (angle of incident acoustical wave relative to orientation of target), shape of target, heterogeneity of target, or acoustical frequency.

Because of the general capability of the DWBA, it has been used to study scatterer types having a wide range of complexity, ranging from the simple smooth homogeneous sphere, for which the DWBA has a deterministic solution, to a rough, heterogeneous deformed cylinder whose randomized boundary and randomized material properties cause the phase shifts of the acoustical waves scattered along the length of the cylinder to vary randomly, as first presented in Stanton et al. (1998). Stanton et al. (1998) showed that randomizing the radii of the infinitesimal cylinder slices and their material properties results in a randomization of the local phase shifts. This statistical formulation relates the randomized phase shifts to variability in scatterer shape, size, and material properties, and the acoustic frequency and angle of incidence. In this way, the parameters of a DWBA model were adjusted to fit scattering data from a live euphausiid over a large range of incidence angles, including those near end-on incidence (Sec. II B of Stanton et al., 1998). In most cases through 2000, applications of the DWBA have involved variations of shape and material properties along the axis of a deformed cylinder. That is, it was incorporated into the DCM (review in Stanton and Chu, 2000). In contrast to these 1-D (line integral) cases, some later applications of the DWBA have involved three dimensions to more accurately account for the scatterer properties (3-D boundary in Lavery et al., 2002; 3-D boundary and 3-D heterogeneities in Jones et al., 2009).

For the uses of the DWBA model in this paper, integration of either the 1-D integral (Stanton et al., 1998; DCM) or 3-D integral (general geometry) (Chu et al., 1993) was performed numerically using the digitized target shapes. Digital resolution and the element length along the axis was 0.1 mm. The DWBA was used to calculate TS dependence on f and θ for the weakly scattering sphere, prolate spheroid, and cylinder.

The PT-DWBA model is a numerical approach for implementation of the DWBA for heterogeneous organisms (Jones et al., 2009). The advantage of this approach is that sound speed and density can vary arbitrarily in three dimensions. For example, material properties can vary due to anatomy of the organism (e.g., different internal organs) and/or irregular external morphology where the acoustic wave passes through part of the animal's body, into the surrounding medium, and back into the body. Thus, the model is conducive to incorporation of high-resolution measurements [e.g., spiral computerized tomography (SCT) scans] of interior and exterior morphology of an organism, and known material properties of the various tissues. The range of acoustic frequencies for which valid predictions can be made using this model is only limited at the high end by the spatial resolution of the morphological measurements.

In the PT-DWBA model the DWBA 3-D volume integral model (Chu et al., 1993; Stanton et al., 1993) is solved numerically for an inhomogeneous medium by integration of the scattering amplitude from each volume element using local material properties to calculate both phase and amplitude of the backscattered signal. The PT-DWBA accounts for the phase change of the acoustic wave by piecewise integration of the phase term along a linear path through the scattering object and surrounding medium.

Three-dimensional, binary matrices were generated in which each matrix element represented a discretization of the scattering volume. Nonzero-valued elements represented the scatterer while zero-valued elements represented the surrounding seawater. The PT-DWBA numerical model was then applied to each simple shape following the model outlined in Jones et al. (2009).

The SDWBA method (Demer and Conti, 2003; Demer and Conti, 2004; Conti and Demer, 2006) is a variant of the statistical representation of the DWBA presented in Secs. II B and II C in Stanton et al. (1998) and summarized above in Sec. III B 4. The two statistical approaches differ in their calculations of random phase shifts. In contrast to the approach in Stanton et al. (1998) where phase shifts are a result of random perturbations in cylinder roughness, material properties, and orientation, the SDWBA creates a random phase shift through a single term that collectively accounts for all sources of phase variability. Both approaches have been shown to improve the prediction of measured scattering levels for angles of incidence well away from normal incidence.

The SDWBA was used to calculate TS values for the weakly scattering sphere, prolate spheroid, and finite cylinder. For each target, the SDWBA results were calculated over all incidence angles for 100 Monte Carlo realizations. The SDWBA is characterized using parameters for standardized length (L₀), number of infinitesimal cylinder slices (N₀), standard deviation of phase variation (φ₀), and frequency (f₀). For this study, values of 38.35 mm, 100, $2/32$⁠, and 120 kHz were used, respectively, except for the finite cylinder, which used N₀ = 50. To predict scattering by idealized (smooth) objects in the absence of variations in target shape and morphology, the value of φ₀ was chosen arbitrarily small to represent negligible roughness.

To maintain a constant ratio of the number of cylinders per wavelength, the cylinder locations for each object were redefined for each frequency. Parameterizations of each object were dependent on the length and the number of cylinders used. The (deterministic) DWBA and SDWBA models are defined by an N + 1 × 3 matrix, xpos, defining the x-, y-, and z-axis coordinates for the centers of the cylinder boundaries, and an N + 1 vector, a, giving the radius at each cylinder boundary. Because the shapes are always assumed to be centered on the z-axis, the z-values in xpos are always set to zero. The parameterizations for each shape are then used to calculate and return the xpos and a arrays for use in the SDWBA.

The physical dimensions of the targets were selected to approximate biological targets such as adult or juvenile stages of fish with gas-filled swimbladders [e.g., Atlantic herring (Clupea harengus)] and zooplankton [e.g., decapod shrimp (Palaeomonetes vulgaris)]. The sphere had a radius of 0.01 m. The prolate spheroid had a minor-axis radius of 0.01 m and a major-axis radius of 0.07 m. The finite cylinder had a radius of 0.01 m and a length of 0.07 m (making the cylinder half the length of the spheroid). The spherical fluid shell had a radius of 0.01 m with a shell thickness of 0.001 m. In addition, depth dependence is ignored (i.e., all targets are at the surface) in these comparisons.

Water temperature and salinity values were chosen to represent conditions in the Gulf of Maine during autumn. Water temperature (6.9 °C) and salinity (34.2 parts-per-thousand) values were representative of the environment north of Georges Bank in 2007. Corresponding water density, ρ_w = 1026.8 kg m⁻³, sound speed, c_w = 1477.4 m s⁻¹ (Table II), at a depth of 0.0 m (i.e., surface adapted), were calculated using Chen and Millero (1977) and Fofonoff and Millard (1983) using the web-based calculators from Chapman (2006) and Tomczak (2000).

Material properties of the targets were selected to approximate those of marine organisms (Table II). For the weakly scattering versions of the sphere, prolate spheroid, and finite cylinder shapes, the density and sound speed were chosen to give g and h values of 1.002. For the weakly scattering spherical shell, the shell density and sound speed were chosen to give g and h values of 1.002 for the shell/water and shell/interior fluid interfaces.

For each target shape and boundary condition, TS was calculated as a function of insonifying frequency and where appropriate, angle of incidence (θ). The frequency range extended from 12 to 400 kHz, which encompasses nearly all bio-acoustical applications at traditional fisheries echosounder frequencies using continuous wave (CW) pulses. Although this frequency range may include resonance scattering for larval fish and other small gas-bearing organisms (e.g., siphonophores), resonance scattering was not included or represented in these comparisons. The frequency response of TS was calculated in 2 kHz increments, where a narrowband CW pulse and implicit time dependence [e^(−iωt)] are assumed. For models that are computationally intensive, TS was computed at f = 12, 18, 38, 70, 120, 200, and 400 kHz. For the prolate spheroid and finite cylinder targets, TS was also computed over a range of angles in the tilt plane from θ = 0° to 90°, where 0° is end-on incidence and 90° is broadside incidence.

TS values were rounded to the nearest 0.1 dB re 1 m² as rounding errors can preclude greater levels of precision. This level of precision is in agreement with practical calibrations of fisheries acoustic instrumentation (Foote, 1983).

The mean, absolute deviation of each non-benchmark TS prediction from its benchmark prediction was used as an index of the match between prediction and the benchmark. This index was computed by taking the mean of the absolute value of the difference in TS between each model prediction and the benchmark at each frequency interval or each angle of incidence.

where N is the total number of frequency or angle values.

Table III and Figs. 1 and 2 present the comparisons between each benchmark and the corresponding model predictions in tabular and graphical format. In addition, the benchmark TS values are given in tabular form in Supplementary Tables I–VI.¹

TS predictions from the BEM were equivalent to the benchmark's TS values from 12 to 400 kHz [Figs. 3(A) and 1(A); Table III]. TS values from the FMM and FEM models were similar to the benchmark, and the KA TS values followed the general trend and amplitude of the benchmark curve, but the values were not equivalent. For the FEM, TS values were nearly identical to the benchmark up to about 80 kHz, and within 0.8 dB up to 200 kHz. For the KA, TS increased from 12 to 100 kHz and then fluctuated about a fairly constant magnitude from 100 to 400 kHz. The BEM, FMM, FEM, and benchmark TS values showed several peaks and nulls across the entire frequency range. Amplitudes of the fluctuations decreased with increasing frequency where they were approximately 4 dB in amplitude at lower frequencies and less than 1 dB at higher frequencies. The KA followed the general trend of the target's backscatter response curve and matched the overall TS magnitude, but had less than half of the peaks and nulls as the benchmark over the frequency range. The KA had a higher maximum TS (−44.0 dB) than the benchmark (−45.4 dB) at about 60 kHz, but had the same overall magnitude from 100 to 400 kHz.

The BEM predictions were equivalent to the benchmark's TS values from 12 to 400 kHz, and the FMM and FEM were nearly equivalent to the benchmark over this frequency range [Figs. 3(B) and 1(B); Table III]. The TS maximum of −42 dB occurred at 12 kHz and decreased monotonically to about −46 dB at 100 kHz, where the TS leveled off at a fairly constant magnitude up to 400 kHz. The TS response curve above 100 kHz was smooth and flat, lacking the peaks and nulls of the rigid sphere.

The BEM predictions were equivalent to the benchmark's TS values from 12 to 400 kHz [Fig. 3(C) and 1(C); Table III]. TS values from the FMM and FEM were nearly identical to that of the benchmark's, but with deviations occurring in the sharp nulls. The general form of the curves was similar to the pressure-release sphere with maximum TS at 12 kHz and a decrease to a constant amplitude from 100 to 400 kHz. There were several sharp peaks and nulls in the backscatter response curves with amplitudes of less than about 2 dB. The KRM TS curve followed the general trend of the target's response curve and matched the overall TS magnitude, with broader peaks and nulls. In addition, the KRM curve showed the opposite trend relative to the other models below about 50 kHz, but did tend toward constant TS at the higher frequencies and a decrease in the amplitudes of the peaks and nulls at higher frequencies.

The BEM, FEM, KRM, DWBA, PT-DWBA, SDWBA, and FMM predictions were all similar to the benchmark's TS values over the 12 to 400 kHz frequency range [Figs. 3(D) and 1(D); Table III]. The curves had broad peaks with narrow, sharp nulls that were 20 to almost 50 dB lower than the peaks. The magnitude of the peaks was fairly constant at −93 dB over the full frequency range. The BEM TS predictions were equivalent to the benchmark up to approximately 120 kHz. At higher frequencies, the mismatch occurred at the nulls where the deviation between the BEM and the benchmark was greater at higher frequencies (approximately 1 dB at 165 kHz and approximately 16 dB at 390 kHz). TS values from the other models were similar, with deviations of usually less than 1 dB, to the benchmark at the peaks, but deviated, sometimes considerably, at the nulls. The modal series–based DCM deviated from the benchmark at the nulls, with two large deviations at 240 kHz and 390 kHz. TS values predicted by the PT-DWBA were equivalent to the benchmark up to about 50 kHz and then deviated at the nulls only. The deviations from the benchmark increased with increasing frequency, with deviations ranging from 1 dB at about 50 kHz to over 20 dB at frequencies over 250 kHz. The SDWBA matched the peaks and nulls well over the entire frequency range. TS values by the KRM generally matched the benchmark at the peaks, with deviations of 1–2 dB, but deviated at the nulls with deviations of 2 to 12 dB.

The fixed-rigid spherical shell is acoustically equivalent to the fixed-rigid solid sphere and therefore a benchmark was not calculated for this case.

The BEM predicted TS values equivalent to the benchmark over the full frequency range of 12 to 400 kHz [Figs. 4(A) and 1(B); Table III]. The BEM predicted monotonically decreasing TS values from 12 to about 100 kHz and then a fairly flat TS response from 100 kHz to 400 kHz, which is similar in form to the pressure-release sphere [Fig. 1(B)]. TS values varied about 4 dB over the frequency range, with the maximum TS (−42.8 dB) at 12 kHz.

The BEM predicted TS values equivalent to the benchmark for the frequency range 12 to 400 kHz [Figs. 4(B) and 1(C); Table III]. The BEM predicted monotonically decreasing TS from 12 to approximately 100 kHz and then relatively constant TS, with the exception of sharp peaks and nulls of 1 to 3 dB in amplitude from 100 to 400 kHz, which is similar to the gas sphere [Fig. 3(C)]. The KRM predicted TS values that followed the general trend of the benchmark TS values above about 50 kHz but with an overall TS level that was 1–2 dB greater than the benchmark. In contrast to the benchmark, the KRM predicted an increase in TS from 12 to about 50 kHz rather than a decrease in TS. The KRM TS values showed several peaks and nulls, but these were smoothly varying as opposed to the sharp spikes in the benchmark.

The KRM, DWBA, and BEM models predicted TS values from 12 to 400 kHz [Figs. 4(C) and 1(D); Table III]. BEM predictions were equivalent to the benchmark's TS values up to about 350 kHz. TS predictions from the PT-DWBA were nearly equivalent to the benchmark from 12 to about 350 kHz, then the DWBA predictions diverged from the benchmark. These models predicted a general decrease in TS over the frequency range with several peaks and nulls. Deviations from the predictions and the benchmark were primarily at the nulls. The KRM predicted similar TS values to the benchmark, but showed an overall increase in TS above 250 kHz compared to the benchmark. In addition, the nulls of the KRM were generally not as deep as those of the benchmark.

Benchmark TS values at broadside incidence for the rigid prolate spheroid were predicted from 12 to 80 kHz [Fig. 5(A) and 1(A); Table III]. The FEM predicted TS values from 12 to 60 kHz, and the BEM, modal series–based DCM, FMM, and KA models predicted TS values over the full frequency range [Fig. 6(A) and 1(A); Table III]. The BEM and modal series–based DCM predictions were equivalent to the benchmark's TS values from 12 to 80 kHz. These TS response curves increased by approximately 7 dB from 12 to 50 kHz then leveled off to slowly undulating TS values. These curves had peaks and nulls with 2–3 dB amplitude at the lower frequencies and less than 1 dB at the higher frequencies. The FMM and FEM TS values followed the same trend and had similar magnitude as the benchmark. The FMM and FEM were similar to each other from 12 to 20 kHz, diverged from each other between 20 and 60 kHz, and then were nearly equivalent from 60 to 80 kHz. The KA had similar TS magnitude as the benchmark, but had only one major peak below 80 kHz. The maximum TS of the KA (−24.6 dB) was greater than the benchmark (−28.6 dB).

Above 80 kHz, TS values predicted by the modal series–based DCM and BEM were nearly identical up to 300 kHz. Above 300 kHz, TS values predicted by the modal series–based DCM and BEM diverged from each other by less than 1 dB. TS values predicted by the KA had similar magnitude to TS values predicted by the modal series–based DCM and BEM, but with fewer peaks and nulls. Peaks and nulls in all three models decreased in amplitude with increasing frequency.

Benchmark TS was predicted from θ = 0° (end-on) to 90° (broadside) from 12 to 80 kHz. The response curves at 38 kHz are shown in Fig. 7(A) and mean deviations in Fig. 2(A). The TS values increased from end-on to broadside incidence with an approximately 40 dB increase and maximum TS of approximately −28 dB at broadside. Amplitudes of the peaks to nulls decreased from nearly 20 dB at end-on to less than 1 dB at broadside. The BEM predictions matched the benchmark at all angles. The modal series–based DCM, KA, FEM, and FMM predictions generally matched the benchmark at angles close to end-on and were nearly equivalent to the benchmark at angles near broadside. The deviation in TS between the DCM and the benchmark increased at angles of incidence near end-on, with the greatest difference of approximately 50 dB at end-on incidence. The FMM had higher TS values than the other models at end-on incidence. The KA followed the pattern of peaks in the benchmark response curve, with general agreement between the peaks of the benchmark and the KA. While KA TS values increased from end-on to broadside incidence, it was generally flat with slight, <1 dB, amplitude fluctuations.

Benchmark TS values from the pressure-release prolate spheroid were predicted from 12 to 80 kHz at broadside incidence [Fig. 5(B) and 1(B); Table III]. FEM predictions ranged from 12 to 60 kHz, and the BEM, FMM, and modal series–based DCM predictions were predicted for the full frequency range [Fig. 6(B) and 1(B); Table III]. The TS response curves were fairly flat with a range in TS values of less than 4 dB from 12 to 80 kHz. BEM predictions were equivalent to the benchmark over 12 to 80 kHz. TS values predicted by the modal series–based DCM were equivalent to the benchmark from 60 to 80 kHz, and only deviated from the benchmark by less than 1 dB from 12 to 60 kHz. TS values from the FMM were less than the benchmark, but only deviated by less than 2 dB of the benchmark's TS values. TS values from the FEM showed two peaks between 12 and 60 kHz, with the FEM predicting the highest TS values.

At frequencies greater than 80 kHz, the BEM and modal series–based DCM models predicted fairly constant TS with a magnitude of about −29 dB. The BEM curve was flat, while the modal series–based DCM curve had several peaks and nulls with amplitudes of about 0.25 dB.

Benchmark TS values for the pressure-release prolate spheroid were predicted from θ = 0° to 90° from 12 to 80 kHz. The response curves at 38 kHz are shown in Fig. 7(B) and mean deviations in Fig. 2(B). BEM predictions were equivalent to the benchmark at all angles of incidence. The benchmark and BEM TS values showed a monotonic increase in TS of nearly 25 dB from end-on to broadside to a maximum of approximately −29 dB. Over the range of angles, TS values from the FMM were greater than the benchmark. The form of the FMM response curve was similar to that of the benchmark, but with greatest deviation from the benchmark at end-on incidence. The FMM matched the benchmark within a few degrees of broadside incidence. The FEM predicted a sharp null at about θ = 10°, but then followed the trend of the other models. In contrast to the other models, the modal series–based DCM had several peaks and nulls. Amplitudes of the fluctuations decreased from end-on to broadside incidence. The modal series–based DCM TS values matched the benchmark at θ > 80°, with the deviation greatest by almost 20 dB at end-on incidence.

No benchmark TS values from the gas-filled prolate spheroid were predicted [Figs. 5(C) and 1(C); Table III]. TS values for the BEM and FEM were predicted to 200 kHz, and the full frequency range for the modal series–based DCM and KRM [Fig. 6(C) and 1(C); Table III]. At frequencies greater than 40 kHz, the BEM, modal series–based DCM, and KRM showed a relatively flat TS response curve, with an overall magnitude of about −29 dB. The KRM response curve had four peaks and nulls of less than approximately 0.5 dB, whereas the DCM response curve had many sharp peaks and nulls of the equivalent magnitude.

TS values predicted from θ = 0° to 90° at 38 kHz are shown in Fig. 7(C) and mean deviations in Fig. 2(C). The KRM and modal series–based DCM predictions were similar to each other with each curve having peaks and nulls in the same angular locations and TS increasing from end-on to broadside incidence. TS values from the modal series–based DCM were consistently less than the KRM, but the difference decreased near broadside to within 2–3 dB at end-on and less than 1 dB at broadside. BEM predictions showed a fairly flat response from θ = 0° to 30°, then a fairly monotonic increase in TS toward broadside incidence. FEM predictions had similar magnitude to the BEM, but with shallower peaks and nulls at angles near end-on.

Benchmark TS values from the weakly scattering prolate spheroid were predicted from 12 to 300 kHz [Fig. 5(D) and 1(D); Table III], and 12 to 400 kHz [Fig. 6(D) and 1(D); Table III]. The benchmark TS values had several peaks and nulls with peak TS values occurring consistently at about −77 dB over the full frequency range. The nulls were narrow and sharp whereas the peaks were broad. In general, the BEM, FEM, modal series–based DCM, KRM, PT-DWBA, DWBA, and SDWBA predictions were similar to the benchmark at the peaks, but deviated, sometimes considerably, at the nulls. TS values by the modal series–based DCM were nearly equivalent to the benchmark over the full frequency range, with only slight deviations of 1 to 3 dB at the nulls. TS values from the BEM were similar to the benchmark, but deviated from the benchmark at the nulls. TS values from the FEM were similar to the benchmark, but also deviated from the benchmark at the peaks as well as the nulls. TS values from the KRM closely matched the benchmark but with deviations of 3 to 20 dB at the nulls. TS values from the PT-DWBA were predicted from 12 to 300 kHz. The TS response curves of the PT-DWBA and DWBA were similar to the benchmark at the peaks, but deviated at the null where deviations generally increased with increasing frequency. In contrast, TS values from the SDWBA matched the benchmark well at the peaks and nulls.

Benchmark TS values were predicted from θ = 0° to 90°. The TS response curves at 38 kHz are shown in Fig. 7(D) and mean deviations in Fig. 2(D). Benchmark TS values increased from end-on to broadside by almost 30 dB to a maximum of approximately −77 dB, and had several peaks and nulls in the response curve. The peaks were broad, whereas the nulls were narrow and sharp. TS values predicted by the PT-DWBA, DWBA, and SDWBA models were very similar to the benchmark TS values over the full angle range, with deviations observed at the nulls. TS values predicted by the BEM were also very similar to the benchmark, with deviations at the nulls. TS values from the modal series–based DCM matched benchmark TS values at angles greater than about 60° and deviated by almost 60 dB at end-on incidence. TS values from the KRM matched benchmark TS values at θ > 60° and deviated from it at end-on by almost 40 dB. Nulls predicted by the KRM did not match the benchmark at θ < 50°.

The benchmark (modal series–based DCM), BEM, and FEM TS response curves for the fixed-rigid finite cylinder were nearly equivalent at frequencies less than about 100 kHz, and equivalent at frequencies greater than approximately 100 kHz [Fig. 8(A) and 1(A)]. TS generally increased from about −37 dB at 12 kHz to −20 dB at 400 kHz. Peaks and nulls in the BEM, FEM, and modal series–based DCM TS response curves had greater amplitudes at lower frequencies. The KA matched the magnitude of the other two models with a monotonic increase in TS over the full frequency range.

TS values at 38 kHz were predicted from θ = 70° to 90° to compare to the benchmark [Fig. 9(A)] and from θ = 0° to 90° [Fig. 10(A)], and mean deviations are shown in Fig. 2(A). TS values from the BEM, FEM, and KA were similar to the benchmark from broadside to about θ = 75°. Between θ = 70° and 75° the angular location of the null differed among the models. Over the full angular range, the predictions were similar to each other, with a slight decrease in TS from θ = 90° to 80°, then a substantial increase in TS approaching broadside. The modal series–based DCM matched the others at θ > 70°, but predicted lower TS values (up to 60 dB lower at end-on incidence) at angles close to end-on incidence. The modal series–based DCM showed a general increase in TS from end-on to broadside incidence. The nulls in the modal series–based DCM response curve were deeper than in the others, especially at θ < 60°.

TS values from the benchmark (modal series–based DCM) and BEM were equivalent for the pressure-release finite cylinder at frequencies above approximately 50 kHz [Fig. 8(B) and 1(B)]. Below 50 kHz, benchmark and BEM TS response curves were within 1 dB of the FEM response curve. All three models predicted monotonically increasing TS over the frequency range, with an increase of about 12 dB to a maximum of approximately −22 dB. The FEM TS response curve showed a similar trend to the BEM and benchmark, but had small, on the order of 1 dB amplitude, undulations.

TS was predicted from θ = 70° to 90° for the benchmark [Fig. 9(B)] and θ = 0° to 90° at 38 kHz [Fig. 10(B)], and mean deviations are shown in Fig. 2(B). TS values were nearly equivalent for the benchmark, BEM, and FEM from about θ = 76° to 90°. The null at approximately 74° was predicted deeper by the benchmark than for the BEM and FEM. Similar to the rigid cylinder, TS values predicted by the BEM and FEM were nearly equivalent over the range of angles. TS values from the BEM and FEM had a slightly decreasing trend from 0° to about 80° with a sharp increase in TS near and at broadside incidence. The modal series–based DCM predicted an increasing trend in TS from end-on to broadside with several deep nulls. The modal series–based DCM predicted lower TS values (up to 40 dB lower) at θ < 60°. At θ > 80° (near broadside) all three models predicted equivalent TS values to a maximum of approximately −31 dB.

The BEM predicted TS values from 12 to 200 kHz, the FEM from 12 to 82 kHz, and the KRM and the benchmark (modal series–based DCM) predicted TS over the full frequency range [Figs. 8(C) and 1(C); Table III]. The benchmark and BEM were equivalent for frequencies greater than about 80 kHz. Below 80 kHz, TS values predicted by the BEM and FEM were slightly greater (less than 1 dB) than the benchmark. The KRM predicted TS values a few decibels greater than the other two models. At frequencies greater than 120 kHz, the KRM and benchmark TS values were within 1 dB of each other. All four models predicted a monotonic increase in TS over the frequency range. The KRM predicted a monotonic increase over the full frequency range, whereas the benchmark TS response curve contained sharp peaks and nulls, of amplitude about 2 dB, at higher frequencies.

TS at 38 kHz was predicted from θ = 70° to 90° for the benchmark [Fig. 9(C)] and from θ = 0° to 90° [Fig. 10(C)], and mean deviations are shown in Fig. 2(C). The BEM, FEM, and KRM were similar to the benchmark from about θ = 75° to 90°. Between θ = 70° to 75° the benchmark predicted a deeper null than for the KRM, BEM, and FEM models. The KRM and modal series–based DCM predicted similar TS values from θ = 30° to 90°. At angles near end-on, the KRM predicted greater TS values than the modal series–based DCM, but both were 10 to 15 dB lower than the BEM and FEM predicted TS values. The KRM and modal series–based DCM predicted generally increasing TS, with several deep (greater than 20 dB) nulls, to a maximum of approximately −31 dB at broadside incidence. The BEM and FEM predicted nearly equivalent TS over the range of angles. The two models predicted slightly decreasing TS from θ = 90° to about 80°, then a sharp increase in TS at broadside incidence. They had shallow (about 5 dB) nulls, which approximately matched the locations of the nulls of the other two models. All four models had nearly equivalent TS values at angles within about 10° of broadside.

The TS values of the weakly scattering finite cylinder were predicted for 12 to 400 kHz by the BEM, 12 to 200 kHz for the FEM, and 12 to 400 kHz for the KRM, PT-DWBA, and benchmark (modal series–based DCM) [Fig. 8(D) and 1(D); Table III]. The BEM, FEM, KRM, benchmark, DWBA, PT-DWBA, and SDWBA models all predicted similar TS values over the frequency range. The TS response curves had several peaks and deep nulls with the nulls being sharp and the peaks broad. The apex of the peaks increased in TS as the frequency increased to a maximum TS of approximately −70 dB. The TS values were similar at the peaks with the greatest deviations at the nulls.

TS values at 38 kHz were predicted from θ = 70° to 90° [Fig. 9(D)] and from θ = 0° to 90° [Fig. 10(D)], and mean deviations are shown in Fig. 2(D). The BEM, FEM, KRM, PT-DWBA, and SDWBA models predicted TS values similar to the benchmark from θ = 70° to 90°. The BEM, FEM, PT-DWBA, and SDWBA models predicted nearly equivalent TS values over the range of angles. TS response curves showed a decreasing trend in TS from end-on to near broadside, and then a sharp increase in TS at broadside. The KRM and modal series–based DCM predictions were similar to the other models from θ = 40° to broadside incidence. The KRM and modal series–based DCM predicted lower TS values, nearly 60 dB, near end-on, but were similar at θ > 40°. All models showed the same deep nulls and broad peaks with similar TS values at θ > 60°.

Few studies have collectively compared predictions from analytical and numerical backscatter models using standard data sets, with exceptions including McClatchie et al. (1996), Foote and Francis (2002), Zampolli et al. (2009), and Macaulay et al. (2013). Most published backscatter model predictions are evaluated by comparison to empirical TS measures (e.g., McClatchie et al., 1996, and references therein), but anatomical differences among animals may increase variability within and among data sets, which may confound results. To illustrate by example, initial comparisons of model predictions by the authors in the current study used anatomical representations of Atlantic cod (Gadus morhua), Atlantic herring, and decapod shrimp as representative gas-filled (i.e., physoclist and physostome) swimbladder-bearing fish and zooplankton. When predicted TS values differed between models, we could not decouple contributions of anatomical scattering structures or model effects from differences in analytical or numerical model structure or their algorithmic implementations. Because of these ambiguities and uncertainties, comparisons based on simple targets were needed to eliminate, or at least minimize, effects of anatomical features on models and predictions. Differences in model predictions could then be attributed solely to the models and their implementation.

While TS values predicted by the BEM matched benchmarks better than the other models, the numerical models (BEM, FEM, and FMM) did not consistently out-perform the approximate analytic models for all target types. For the rigid, pressure release, and gas spheres, and the rigid and pressure release prolate spheroid targets, the BEM was within 0.1 dB average deviation of the benchmarks across the modeled frequency and angle of incidence ranges. Predictions from analytic approximations, such as the KA, modal series–based DCM, and KRM compared favorably with those of the benchmark for spherical targets (KA only) and broadside incidences of the prolate spheroid, where most were within 1 dB average deviation over the 12 to 400 kHz frequency range. Even though overall trends and magnitudes of the approximations matched the benchmarks, TS was not well predicted at the sharp peaks and nulls that are a feature of many of the benchmark backscatter response curves. The “exact” or “highly” accurate models (Table I) had more structure (i.e., sharp peaks and nulls), where each feature ranged over a few kHz, and TS amplitude variations did not follow the overall TS trend. In contrast, the KA-based models predicted gradually undulating variations where the frequency of the variation ranged 10 kHz or more, and amplitudes of the variations tended to decrease with increasing frequency.

In general, predicted TS values closely matched those of the benchmarks at angles near broadside incidence, but did not match well at off-broadside incidences. For the prolate spheroid, average deviations ranged from 1 to over 9 dB when off-broadside angles were included. TS predictions from the modal series–based DCM and KRM models deviated by 20 dB or more at end-on incidence, but were nearly equivalent to the benchmark values at or near broadside incidence. This is expected for those two models as they are inherently inaccurate at angles well off normal incidence. For the cylindrical targets, the modal series–based DCM was used as the benchmark for broadside incidence and for angles from 20° off-broadside to broadside incidence. The choice of 20° was based on investigations of rigid cylinders by Partridge and Smith (1995) for ka ranges of 10 to 40 and aspect ratios between 2.5 and 24.0. Ye (1997) demonstrated that end effects become noticeable when the angle of incidence was greater than 40° at ka = 40 and an aspect ratio of 2.55. From results presented here, it appears that the modal series–based DCM is a valid benchmark for angles up to approximately 15° for rigid, pressure-release, and gas-filled finite cylinders, and up to approximately 20° (from broadside) for a weakly scattering finite cylinder.

Maximum TS values and overall trends of the weakly scattering targets matched those of the benchmark, but locations and amplitudes of the predicted nulls were often offset from and did not match the benchmark. Both of these effects caused large, often greater than 10 dB, differences between predicted and benchmark TS values. The locations of the nulls are dependent on the material properties and orientations of the targets (Stanton and Chu, 2000). These results also show that null locations are dependent on the numerical or analytical models used to predict TS values. The DWBA, PT-DWBA, and SDWBA models predicted similar TS values over all frequency and angular ranges, although the TS values were not equivalent. Differences in predicted TS values suggest differences in the implementation of the algorithms since they are mathematically equivalent when parameterized for these particular targets. That is, the targets used in this exercise do not incorporate any additional anatomical complexities of the marine organisms that the three models were developed to address, such as boundary roughness and heterogeneities. Differences between predictions for the more complex scatterers would be expected due to the different approximations made in each of the models, as discussed above.

Most of the models included in the comparisons required a digital representation of target shape, but the level of detail (i.e., resolution) required to predict TS differed among the models. The numerical models are sensitive to the resolution of the digital representation, especially for convergence. The BEM, FEM, and FMM models require discretization of the model domain into elements that must be sufficiently small relative to the acoustic wavelength (∼λ/10) to obtain an accurate or even an adequate solution (Stanton and Chu, 2000; Foote and Francis, 2002). As frequency increases, wavelength decreases, and element size is reduced, computational and computer memory requirements increase, especially for 3-D representations used in this exercise. Predictions from the FEM and FMM models were within 1 dB average deviation over the 12 to 400 kHz frequency range of the benchmark for spherical targets and broadside incidences of the prolate spheroid. For the prolate spheroid and fluid shell targets in FEM modeling and for most targets at the higher frequencies it was not possible to maintain the necessary target discretization resolution, and model accuracy was correspondingly degraded, and/or the maximum frequency was restricted to lower values.

A trade-off exists between computational requirements and prediction accuracy. The BEM is a superior model for predicting TS values for these benchmark shapes, but is computationally intensive. Approximations such as the KA and KRM models do not require extensive computational resources, but accuracy can be reduced, especially with increasing off-broadside incident angles (Macaulay et al., 2013). It should also be noted that the resolution of the digital representation of the target can affect the amplitude and variability of predicted TS (Jech and Horne, 1998; Stanton and Chu, 2000). This was especially an issue with use of the KRM that had a fixed resolution (as per its operational use) which was too coarse at the higher frequencies. Although the trends of predicted TS from the approximations were frequently within 1 dB of the benchmarks, they did not always match the fine-scale structure present in benchmark backscattering response curves.

TS predictions by analytical and numerical models and digital representations of organism anatomy are becoming essential tools in fisheries acoustics research. They are used to predict TS values over a range of conditions that cannot be practically sampled in the field (Reeder et al., 2004), to estimate TS variability across life history stages (Horne, 2008), for aggregations of animals, and to determine the relative importance of physical and biological factors influencing TS. Selection of a single model for all types of organisms and environmental conditions is unrealistic given the limitations of each model for scattering types, computational intensity, and required resolution of the digital representation. For the targets presented here, the numerical models worked well at the lower ka ranges (e.g., < 5), but as the ka range increased and the target shape became more complex, the computational load increased such that predictions required hours to days to compute. Analytical models and approximations were accurate within the limitations of specific target types, target orientations, and ka ranges in this study, but these ranges can be restrictive when applied to animals.

At present, a suite of models is required to predict acoustic backscatter from aquatic organisms with complex shapes, morphologies, and behaviors. In this exercise, targets were idealized using simple shapes and morphologies; such as weakly scattering targets representing zooplankton and the gas-filled targets representing gas-filled bubbles such as swimbladders. All of the numerical and analytical models for weakly scattering spheres generally matched the benchmarks, which was also the case for weakly scattering spheroidal and cylindrical targets at or near broadside incidence. This suggests that approximate analytical models may be preferred over numerical models for these targets as they are computationally less intensive and provide accurate TS predictions for, at least, these idealized objects. For TS predictions at incident angles well off broadside, numerical models may be preferred even though they can be computationally intensive. For weakly scattering targets, the (analytical) DWBA-based models performed well for these targets. For gas-filled targets, numerical models were superior to approximate analytical models for exact matching of TS values over the range of frequencies and angles examined. In most cases, with the exception of off-broadside incidence, approximate analytical models predicted overall TS trends, and may be sufficiently accurate for density and abundance estimates, especially when these estimates are for groups of idealized objects with varying shapes, sizes, and material properties. Approximate analytical models for density and abundance estimates can be computed quickly, do not require specialized software, and often do not require the high morphometric resolution needed by numerical models. Operationally this means that approximate analytical models can be quickly generated in response to changing biological or environmental conditions during, for example, a research survey. But with sufficient time, numerical models provide accurate and precise TS predictions, at least for these idealized objects, that can be applied after data have been collected. The next logical step in the evaluation of predicted backscatter accuracy and precision is to compare TS predictions of real fish and zooplankton using a common set of morphological data and ambient noise conditions, including measures of anatomical variability.

This work was supported by the NOAA Fisheries Advanced Sampling Technologies Working Group, the Office of Naval Research, and the National Oceanic Partnership Program. Josiah S. Renfree, SWFSC, is thanked for computing the SDWBA predictions. K. Foote is thanked for comments and discussions. Two anonymous reviewers are thanked for constructive comments. C. S. Clay and D. Van Holliday contributed not only to the early stages of this work, but also throughout their careers to establishing the foundation that produced this work.

The coefficients in Eqs. (8) and (9) are expressed in terms of the Spherical Bessel functions of the first [j_n(x)], second [y_n(x)], and third [$hn1(x)$⁠, or Hankel function].

Notations are a and b are radii of the outer and inner spheres, respectively, and $gij=ρi/ρj$ and $hij=ci/cj$⁠, where ρ_k and c_k are density of and sound speed in the kth medium. Subscripts 1, 2, and 3 represent the surrounding water, outer sphere, and inner sphere, respectively.

1. Coefficients in Eq. (9):

   $b1=jn(k1a),$

   $b2=g21h21jn′(k1a),$

   $a11=−hn(k1a),$

   $a21=−g21h21hn′(k1a),$

   $a31=0,$

   $a12=jn(k2a),$

   $a22=jn′(k2a),$

   $a32=jn(k2b)jn′(k3b)−g32h32jn′(k2b)jn(k3b),$

   $a13=yn(k2a),$

   $a23=yn′(k2a),$

   $a33=yn(k2b)jn′(k3b)−g32h32yn′(k2b)jn(k3b).$
2. Coefficients in Eq. (10): Coefficients b₁, b₂, a₁₁, and a₂₁ are given above.

   $d1=jn(ksa)yn(k2b)−jn(k2b)yn(k2a),$

   $d2=jn′(ksa)yn(k2b)−jn(k2b)yn′(k2a).$