Estimating target strength and physical characteristics of gas-bearing mesopelagic fish from wideband in situ echoes using a viscous-elastic scattering model

Target strength (TS) from gas-bearing organisms with gas inclusions is very distinct from organisms lacking them (Stanton et al., 2010). In the mesopelagic zone (200–1000 m) some species of fish [e.g., Marshall (1960), Butler and Pearcy (1972), Davison (2011), and Scoulding et al. (2015)] and physonect siphonophores [e.g., Barham (1963), Kloser et al. (2016), and Proud et al. (2019)] have gas inclusions and are therefore strong scatterers of sound (Foote, 1980; Leighton, 1994).

Mikronekton (2–20 cm in size) organisms inhabiting the mesopelagic zone are known to be numerous, and the mesopelagic fish component of this is estimated to comprise a high fraction of the estimated total global fish biomass (Irigoien et al., 2014). Mesopelagic fish are a potential source of marine fat and protein [e.g., Gjøsæter and Kawaguchi (1980)] and play key roles in ocean ecosystems by transferring energy from lower to higher trophic levels (Beamish et al., 1999) and by contributing to active carbon transport from surface to deep waters through diel vertical migration (Gjøsæter and Kawaguchi, 1980; Irigoien et al., 2014; Proud et al., 2019). Their actual ecological importance is somewhat uncertain, partially due to high uncertainties in biomass estimates, with current global estimates spanning about one order of magnitude (Gjøsæter and Kawaguchi, 1980; Irigoien et al., 2014; Proud et al., 2019). A key in fisheries acoustics is to understand the backscattering of organisms for biomass estimation and/or other quantitative studies (Ona, 1999). Here, we used a two-layer viscous-elastic spherical gas backscattering model to estimate physical parameters in gas-bearing fish by using in situ broadband backscattering measurements.

Observation and quantification of deep-living marine organisms are often attempted using acoustic backscattering measurements (Simmonds and MacLennan, 2005; Irigoien et al., 2014; Proud et al., 2019), which are an efficient high-resolution complement to other techniques such as optical observations and trawl catches. Acoustic sampling is not only more efficient but also a non-extractive method for monitoring of marine organisms. Echograms, which are a visualisation of the acoustic data, provide information about location and signal strength of targets in the water column, but identification of those acoustic targets to species or species groups remains imprecise (Simmonds and MacLennan, 2005). To convert the measured backscattered acoustic signal from marine organisms to biological information (such as size, weight, and species identification), it is necessary to understand how the organisms scatter the incident acoustic waves (Horne, 2000; Reeder et al., 2004), particularly the proportion of acoustic energy a single individual scatters back to the receiver [the backscattered target strength (TS)].

Gas bubbles are strong scatterers of sound because of their large acoustic impedance contrast compared to water (Leighton, 1994). A gas-filled inclusion (e.g., swimbladder or pneumatophore) is the main reflector of sound for species that possess one. For a fish with swimbladder, more than 90% of total backscattered energy is caused by the swimbladder (Foote, 1980). As gas-bearing organisms are strong acoustic reflectors, they will, if misinterpreted as non gas-bearing, result in significantly biased biomass estimates (McClatchie and Coombs, 2005). The target strength frequency response from a swimbladder-bearing species of fish will be very distinct from that of a species lacking a swimbladder, enabling the distinction and categorization of these two different groups (Stanton et al., 2010). The scattering of acoustic waves by bubbles are more pronounced at resonance and are about three orders of magnitude larger than the geometrical scattering (i.e., where acoustic wavelength are much smaller than bubble size) (Leighton, 1994). Many mesopelagic fish species have a gas-filled swimbladder at least in the early stages of their life (Butler and Pearcy, 1972; Yasuma et al., 2003; Yasuma et al., 2010). Also, some species show ontogenetic and latitudinal differences in swimbladder sizes, and fat-invested swimbladders have been observed mainly in larger individuals of a species (Marshall, 1960; Butler and Pearcy, 1972; Davison, 2011; Scoulding et al., 2015; Dornan et al., 2019).

Previous studies suggest that physonect siphonophores contribute to backscattering from within the mesopelagic zone in some areas [e.g., Barham (1963), Kloser et al. (2016), and Proud et al. (2019)]. Their gas-inclusion can resonate at depth (Kloser et al., 2016; Knutsen et al., 2018; Proud et al., 2019) at similar frequencies as that for mesopelagic gas-bearing fish (Kloser et al., 2016). The presence of siphonophores has been documented optically [e.g., Kloser et al. (2016) and Knutsen et al. (2018)] and in nets (Greene et al., 1998; Knutsen et al., 2018) and if present in large numbers, this could lead to the overestimation of mesopelagic fish biomass if only fish were assumed present (Proud et al., 2019).

TS depends on an organism's material properties (density and sound speeds), shape, size, and orientation in relation to the incoming sound pulse (Faran, 1951; Hickling, 1962; Stanton et al., 1998). Backscattering modeling enables the study of these factors separately (Hazen and Horne, 2003). Therefore, target strength modeling of individual organisms provides information that can be used to fill the gap between measured backscatter and biology and can reduce uncertainties in abundance estimation. Furthermore, theoretical models can enhance the ability to size, recognize, identify, and discriminate acoustically measured targets (Jech and Horne, 2002).

There are several theoretical models for bubble and swimbladder acoustic backscattering [e.g., Anderson (1950), Love (1978), Feuillade and Nero (1998), and Ainslie and Leighton (2009, 2011)]. Because of the complex shape and inhomogeneous tissue of marine fish and their swimbladders, the backscattering is inevitably approximated by theoretical models of simple geometric shapes (Anderson, 1950; Foote, 1980) such as a spherical fluid-filled bubble (Anderson, 1950), gas-filled cylinder (Clay, 1991), or prolate spheroid (Ye, 1997; González et al., 2016). The spherical shape is the simplest and most widely used to model swimbladder backscattering and has provided significant insight into aquatic organisms (Medwin, 2005). However, to better resemble the structure of a swimbladder, Love (1978) added a viscous-elastic shell to the spherical bubble and used only the monopole backscattering mode. Other models included higher modes to model backscattering from a spherical bubble with elastic (Goodman and Stern, 1962) and viscous (Anson and Chivers, 1993) shells. Feuillade and Nero (1998) modeled resonance scattering of a swimbladder by a gas-filled sphere enclosed by an elastic layer to represent the swimbladder wall, which was surrounded by a viscous layer representing fish flesh. The model by Feuillade and Nero (1998) includes the rigidity and damping effects of the swimbladder wall and fish flesh, respectively, and can be used to provide information on the acoustic and mechanical properties of live tissues, which are scarce especially for mesopelagic fish.

Broadband acoustic backscatter techniques have several advantages over conventional multifrequency narrowband methods (Stanton et al., 2003). Multifrequency backscatter techniques can be used to distinguish between major scattering groups based on their relative frequency response measured at discrete frequencies [e.g., Korneliussen and Ona (2002) and Korneliussen et al. (2016)]. The discrete narrowband frequency response information is usually not sufficiently detailed to separate acoustically similar species or different size groups of a single species (De Robertis et al., 2010). In comparison, broadband acoustic backscatter can provide frequency response over broad frequency intervals (Horne, 2000) that can potentially enhance acoustic identification and provide information on target properties such as morphology or size (Reeder et al., 2004; Kubilius et al., 2020), and has led to improved species discrimination (Stanton et al., 2010; Verma et al., 2017; Bassett et al., 2018). In addition, broadband acoustic backscatter signals can provide significantly increased range resolution compared to narrowband systems through matched filtering (Lavery et al., 2010; Stanton et al., 2010), thereby enabling TS measurements of single organisms in denser aggregations. Previous studies using broadband acoustic backscatter have tended to focus on frequency responses of volume backscattering strengths (S_v), i.e., focusing on aggregations rather than single individuals [e.g., Bassett et al. (2018) and Benoit-Bird and Waluk (2020)]. For mixed assemblages such approaches are anticipated to be inaccurate with the frequency response of weaker scatterers likely to be masked by stronger scatterers.

Although widely used and highly applicable for studying organisms in the upper water column, one limitation of echosounders affixed to the hull of vessels is the limited working range of higher frequencies due to acoustic absorption by water (Francois and Garrison, 1982a,b), which limits vessel-mounted acoustics to about 100 kHz or lower when studying the deeper parts of the mesopelagic zone. Furthermore, the longer the range, the higher the probability of getting more than one organism within the acoustic beam and restricting measurements from individual targets. By lowering acoustic instruments to the depth of interest these issues can be overcome and allows for the use of broadband acoustic signatures for identification of organisms at mesopelagic depths (e.g., Verma et al., 2017). This approach can provide information on both taxonomic composition and densities inside deep scattering layers (DSLs), where densities are poorly known. It is, however, important to note that ground-truthing of acoustic targets remains essential and therefore trawls and/or optical methods are still needed [e.g., Kloser et al. (2016)], although both of these methods have their own problems and inherent biases.

In situ frequency responses combined with a backscattering model have the potential to provide characteristic features of the target and result in a tuned model of mesopelagic gas-bearing organisms. In this paper we use a towed platform equipped with acoustic transducers, deployed within the mesopelagic zone, to overcome the limited range of higher frequencies at depth. We used this platform to collect wideband acoustic measurements from mesopelagic organisms. Based on biological sampling using trawls, photographic images obtained from a stereo camera system (Deep Vision) (Rosen and Holst, 2013) attached to the aft of the macroplankton trawl, and a video plankton recorder (VPR) (Davis et al., 1992) attached to the front of the towed platform, we concluded that the observed gas-bearing targets are most likely to be fish. Furthermore, we used a two-layer mathematical/physical backscattering model that incorporates higher modes of scattering to fit the measured acoustic data and estimated swimbladder sizes. The potential for estimating swimbladder wall thickness, its shear modulus, and flesh viscosity were also investigated.

Data used in this paper were collected during a research cruise in the eastern part of the Mid-Atlantic Ocean from Cape Verde to the southern part of France (17°N 25°W to 48°N 8°W) (Fig. 1) on board R/V Kronprins Haakon (Norwegian Institute of Marine Research, IMR) 2nd to 22nd of May 2019. The objective was to advance our understanding of the mesopelagic ecosystem along latitudinal and longitudinal gradients in the study area.

A towed vehicle (MESSOR) (Knutsen et al., 2013) equipped with a four channel echosounder (Simrad EK80 WBT Tubes operating at 38, 70, 120, and 200 kHz) was used to collect acoustic data from within the mesopelagic zone at 15 stations (Fig. 1). The transducers were mounted on the bottom plate of MESSOR, facing downwards. Three of the depth-rated transducers (70, 120, and 200 kHz) were operated in FM (frequency modulated) mode, whereas the 38 kHz transducer was operated in CW (continuous wave) mode. Interference and crosstalk between channels was minimized by reducing the transmit power of the lower frequency transducers [see Table I and Khodabandeloo et al. (2021) for details]. The listening range was reduced to 60 m as the reduced power also reduced the signal-to-noise ratio. The ping interval was varied as needed to reduce interference from double (surface and bottom) returns on the 38 kHz channel and were kept as small as possible, generally ranging from ∼250 to 350 ms.

MESSOR was towed behind the ship in oblique hauls from 0 to 1000 m depth for 4 h at a speed of ∼2 m s⁻¹. Acoustic data collected with 38 kHz narrowband and 50–80 kHz, 93–155 kHz, and 160–260 kHz broadband (see Table I for data collection settings) were processed to yield in situ measurements of the target strength frequency response from mesopelagic organisms. MESSOR was further equipped with a conductivity, temperature, and depth (CTD) profiler (Seabird SBE 49 FastCAT) that was operated throughout the deployments. The CTD data were used in Eqs. (A13)–(A15) to estimate the densities and sound speed of the surrounding seawater as a function of MESSOR depth. This is explained in Sec. II B 2 b.

Calibration of the echosounder system was conducted using standard methods (Demer et al., 2015) at the surface (∼2 m depth), at a station off Cape Verde on May 2nd 2019. A CTD profile (SeaBird 911plus) was taken prior to calibration of the echosounder system to estimate seawater sound speed. A tungsten carbide (with 6% cobalt binder) 38.1 mm diameter sphere was used for calibration of all four transducers using the same settings as during data collection (Table I).

The acoustic data used in the present study are derived from a MESSOR profile taken from 32.54°N 12.11°W to 33.08°N 12.30°W (Fig. 1). We manually identified single targets at depths ranging from ∼300 to ∼900 m in order to retrieve target strength frequency responses from mesopelagic organisms.

The acoustic data from MESSOR was postprocessed in the lsss computer program (Large Scale Survey System) (Korneliussen et al., 2016). We manually chose individual targets with the criterion that each target had been observed on all four frequency bands. Target strength as a function of frequency was obtained from a fast Fourier transform (FFT) of the pulse compressed echoes using an FFT window length of 0.3 m centred on the target (Table II). Selecting a long window (vertical extent) around a target includes more backscatter information and hence provides higher resolution and more complete frequency response of the target. However, especially in dense scattering layers, long windows around adjacent targets can overlap, which distorts the frequency response due to interference between the backscattered signals from multiple targets (Stanton et al., 1996; Reeder et al., 2004). For example, such interference lead to regularly spaced nulls in a frequency response (see, e.g., “target C” in Fig. 2). In summary, the proper length of FFT window is a trade-off between separating individual targets and gaining sufficient frequency response information from single targets.

If a given target met the desired parameters (Table II), the TS detector accepted it as a target and provided a corresponding TS frequency response. Only targets within the nominal 3° beam angle were included for further analysis.

The correct measurement of TS relies on the selection of echoes from a single target and avoiding multiple echoes from adjacent targets. Hence, to ensure each manually selected target had a high probability of being from a single organism, we used both frequency response and target location information within the acoustic beam (see Fig. 2). As examples of accepted single targets, see targets A and B. These two targets have straight consistent direction (middle panels, Fig. 2) from ping to ping as well as the consistency between TS measurements (lower panels, Fig. 2). Targets C–E are rejected as single targets because target locations (middle panel), null patterns in the frequency responses and inconsistency of TS measurements (lower panel) indicate that there is more than one target inside the beam.

The frequency response data from selected targets were then exported from lsss for further analysis.

The in situ broadband frequency response data were used to fit a parameterized scattering model, which will be described in Sec. II B.

To estimate backscattering from a swimbladder, the mathematical/physical model introduced by Feuillade and Nero (1998) is used. The model (Fig. 3) assumes a spherical swimbladder and includes swimbladder wall tissue using an elastic shell and also the damping effects of fish flesh on the scattered acoustic signal via a surrounding viscous layer.

Feuillade and Nero (1998) expressed the wavefield potentials within the target in terms of spherical Bessel and Neumann functions. Avoiding instability of this formulation at high frequencies (i.e., well above main resonance) for the higher modes (m > 0) is complicated, but these modes are particularly important for TS modeling in the high frequency region, and should therefore be included. A more straightforward method is to use spherical Hankel functions instead of Neuman functions (Anson and Chivers, 1993) and this is further elaborated in Sec. II B 1. We present the model derivation here to provide a self-contained reference for discussion about effects of the input parameters.

The velocity potential of incident plane wave field with unit amplitude is given by

where $ω$ is the angular frequency, $k1(=ω/c1)$ is the wave number in the water, $θ$ is the planar wave incident angle, $Pm$ is the Legendre polynomial of order $m$⁠, and $jm$ is the sphercial Bessel function of the first kind, order $m$⁠.

The compressional wavefield reflected in the water is given by

where $hm$ is the spherical Hankel function of the first kind and order $m$⁠. $Am(1)$ is the coefficient to be determined. The velocity potential of compressional and shear field in the viscous layer are given, respectively, as

The compressional and shear wave numbers in the viscous layer are, respectively,

In the above equations, $μ2$ is the coefficient of shear viscosity and $ζ2$ is the viscosity parameter defined as (Love, 1978; Feuillade and Nero, 1998)

where $η2$ is the coefficient of bulk viscosity in the viscous layer (i.e., medium 2).

The velocity potential of compressional and shear field in the elastic shell are written, respectively, as

The velocity potential of compressional and shear wave numbers in the elastic layer are obtained via

where $λ3$ and $μ3$ are Lamé constants in the elastic medium. For Lamé constants corresponding with soft to hard rubber, the scattering is only sensitive to the shear modulus (⁠$μ3$⁠). Therefore, effects from Lamé parameter $λ3$ are ignored (Feuillade and Nero, 1998) and assumed to be constant.

Finally, the velocity potential of the compressional wave field within the gas inside the swimbladder is given by

where $k4 (=ω/c4)$ is the wave number in the gas and $c4$ is the sound velocity in the swimbladder gas. To estimate the TS around the resonance, only the first term (m = 0), which corresponds to the monopole (i.e., pulsation due to volume change without deviation from spherical shape), is sufficient. But since we want to estimate TS at frequencies well beyond the resonance, more terms need to be included.

The appropriate boundary conditions at the three interfaces between media are

1. Continuity of normal velocity between medium 1 and 2 at $r=R2,$

2. Continuity of normal stress between medium 1 and 2 at $r=R2,$

3. Tangential stress equal to zero at $r=R2$⁠,

4. Continuity of normal velocity between medium 2 and 3 at $r=R3,$

5. Continuity of normal stress between medium 2 and 3 at $r=R3,$

6. Continuity of tangential stress between medium 2 and 3 at $r=R3,$

7. Continuity of tangential velocity between medium 2 and 3 at $r=R3,$

8. Continuity of normal velocity between medium 3 and 4 at $r=R4,$

9. Continuity of normal stress between medium 3 and 4 at $r=R4,$

10. Tangential stress equal to zero at $r=R4$⁠.

From the above list, items 1, 5, 7, and 10 include normal/tangential velocity and stress boundary conditions and are given in the  Appendix. The rest are similar to one of these items and can be expressed accordingly. Applying the abovementioned ten boundary conditions to find ten unknown coefficients $Am1, Am2$⁠, $Bm2$⁠, $Cm2, …$⁠, $Am4$ in Eqs. (2), (3), (7), and (10), a system of equations are written in matrix form as

Note that for the first mode of backscattering (i.e., M = 0) there are six unknowns instead of ten and therefore six proper boundary conditions (in this case boundary conditions 1, 2, 4, 5, 8, and 9) should be used [more details in Feuillade and Nero (1998)].

The pressure in terms of the velocity potential is expressed as (Pierce, 1989)

and the backscattering form function is defined as (MacLennan, 1981; Reeder et al., 2004; Jech et al., 2015)

where $pinc$ and $pscat$ are the amplitude of the incident and scattered pressure at distance $r$ from the target.

Target strength is then given by

Since the incident pressure is known [Eqs. (1) and (12)], the backscattered pressure is required to estimate the TS [Eqs. (13) and (14)]. To estimate the backscattered pressure [Eqs. (2) and (12), $θ=π$⁠)], $Am(1)$ is required, which can be determined by solving the system of equations [i.e., Eq. (11)].

To solve the system of equations given in Eq. (11) care must be taken to avoid instabilities. The instability is caused by the imaginary part of wave numbers which occurs for the shear wave numbers in the viscous medium (Anson and Chivers, 1993). As a result, some of the elements of matrix $H$ are extremely large and therefore cause an ill-posed problem. Writing the wavefields in terms of spherical Hankel - instead of Neuman - functions, the coefficients of $Cm2$⁠, $Dm2, Cm3$⁠, and $Dm3$ in the $H$ matrix are expressed only in terms of spherical Hankel functions. Therefore, it is possible to factor out the exponential terms using the recursive formulas for the spherical Hankel functions of first kind as (Arfken and Weber, 2005)

For the complex values of the wavenumber (the argument) the spherical Hankel function is expressed as

where the imaginary part can be factored out.

The model has a total of twelve parameters (Table III) that are needed to calculate a backscattering frequency response.

The gas inside the swimbladder of mesopelagic fish is mainly oxygen (Ross, 1976; Wittenberg et al., 1980; Priede, 2017). To estimate gas density at the pressure experienced in situ, both Boyle's model [Eq. (A12)] and Van der Waals' model [Eq. (A13)] were used (Priede, 2018).

Calculated sound speed [using Eq. (A14)] in oxygen at atmospheric pressure for three temperatures (T = 1 °C, 13 °C, 30 °C) are shown (solid square in Fig. 4, right panel). Experiments by Van Itterbeek and Zink (1958) show that the sound speed increases slightly with increasing pressure (Fig. 4, right panel). To estimate the sound speed inside the swimbladder within a mesopelagic layer at a given temperature, we have fitted a line to each of the experimental data above 20 bar for two measurements at temperatures 1 °C and 30 °C (Fig. 4, right panel). The sound speed for a desired temperature (1 < T < 30) and pressure (>20 bar) can be found using weighted averaging relative to the temperatures 1 °C and 30 °C. For example, for T = 13 °C the sound speed is shown as function of pressure between 20 and 100 bar (Fig. 4 right panel).

From the CTD mounted on MESSOR we obtained in situ salinity and temperature measurements from surface to 1000 m depth. Measurements between 200 and 800 m were used to estimate density of seawater using Massel (2015), Appendix A. The sound speed of seawater was estimated [see Eq. (A15)] also using in situ measured temperature, salinity as a function of depth.

If the swimbladder makes the fish neutrally buoyant, then

where $ρg$ is the density of gas inside the swimbladder, $ρw$ is the water density, $ρf$ is the fish flesh density, and $Vf$ is fish volume.

We assume that the fish are neutrally buoyant. Neutral buoyancy significantly reduces the required energy for swimming (Priede, 2017). Furthermore, we assume that all the fish flesh is concentrated around the swimbladder as a sphere with equivalent spherical radius (henceforth “radius”) $R2$⁠. Equation (17) is then expressed as

There is little information available about the elastic properties of fish swimbladder tissue, particularly for mesopelagic fishes. The shear modulus of rubber was used to model the cod swimbladder by Feuillade and Nero (1998). They used values from 0.3 to 2 MPa, corresponding to the shear modulus of soft to hard rubber. Sand and Hawkins (1973) estimated the shear modulus of cod swimbladder to be 0.17 MPa but suggested that it can increase to 1–20 MPa at depth. In a recent paper, Fine et al. (2016) measured Young's modulus of oyster toadfish (Opsanus tau) swimbladder wall to be between 0.5 and 3 MPa. The measured values vary for anterior and posterior and different parts of the swimbladder (dorsal, side, and ventral) and furthermore depends on direction (normal or parallel to long axis of the swimbladder). Using Young's modulus, the shear modulus is estimated to be around 0.17 to 1 MPa, assuming the Poisson's ratio of ∼0.5 [0.4999 for rubber, Engineering ToolBox (2008)] for the swimbladder tissue. Even though shear modulus is, in general, a frequency dependent parameter it can be assumed to be constant above 1 kHz [Ye and Farmer (1994) and Fig. 13 in Carstensen and Parker (2014)].

The thickness of swimbladder walls have been measured in some mesopelagic fish species and found to vary from 10 to 300 μm (Marshall, 1960) (Table V). The values used in our model were within this range.

Backscattering from a gas-filled sphere over a frequency of 1–260 kHz was estimated using the finite element method (FEM) to benchmark other solutions. A three-dimensional FE model (adapted from acoustic scattering off an ellipsoid, comsol multiphysics 5.4, Acoustics Module model Library) was implemented using the comsol multiphysics^® v.5.4 computer program. The inhomogeneous Helmholtz equation was solved in the frequency domain, where the physical domain was discretized into small elements such that there were at least 18 elements per wavelength. A perfectly matched layer (Berenger, 1994) was used around the computational domain (water) to attenuate the waves reaching the computational domain boundary in order to minimize the effects of using a finite computational domain. FEM is computationally demanding, especially at higher frequencies (Jech et al., 2015), and at the frequency regions where there were peaks and valleys in the backscattering curves, a finer frequency resolution was used. It was observed that for correct backscattering estimation, the water layer surrounding the target should resolve to be at least a quarter of wavelength thick. For example, for a frequency of 15 kHz, the wavelength in the water is 0.1 m and hence, the water thickness should be at least 0.025 m. To avoid an unnecessarily large model domain, and consequently large computational effort, the water domain was extended to a quarter of wavelength and hence varied at each frequency. The perfectly matched layer thickness was set to a thickness of one-eighth of the wavelength. The far-field backscattered pressure was calculated by solving the Kirchhoff-Helmholtz integral equation using the appropriate Green's function at a sufficiently large distance (⁠$R→∞$⁠), ignoring the rapidly oscillating phase factor (comsol multiphysics 5.4, 2018):

where $S$ is a closed surface (here surface of the water domain) and $n$ is a normal unit vector into the domain enclosed by the surface $S$⁠. The incident pressure amplitude was set to 1 Pa and TS was calculated using Eq. (14).

By making the viscosity and elastic properties of the swimbladder wall very small, the viscous-elastic model presented in this paper should resemble a gas sphere in a fluid. Therefore, the viscous-elastic model can be benchmarked against other existing models for gas bubbles [e.g., modal solution by Anderson (1950)]. In this regard, the viscous-elastic backscattering model was validated and the minimum number of required modes (i.e., terms in the summations) for the target strength estimation determined. Backscattering was calculated using one, two and three terms (Fig. 5 comparing M = 0, 1, and 2). Note that M = 2 means summation of terms with m = 0, 1, and 2 [Eqs. (1), (2), (3), (7), and (10)]. The resulting backscattering from the soft viscous-elastic model including one, two, and three modes were compared to the backscattering from the Anderson modal solution and finite element (FE) method (Fig. 5). In addition, Love's model (Love, 1978), which is widely used in fisheries acoustics for backscattering estimation from a gas sphere, is plotted for the case of zero damping (Fig. 5).

The viscous-elastic model is in good agreement in the main resonance region when including just the first mode (zero order mode, m = 0) which is equivalent to monopole backscattering (i.e., M = 0) (Fig. 5). At higher frequencies, the backscattering is overestimated [e.g., Fig. 5(C) and lower left panel] and the one dip is not present in the other models [e.g., Fig. 5(D)]. When two modes are included (i.e., M = 1), the peak at around 182.6 kHz is missed and there is a slight difference for the dip appearing at around 245.3 kHz compared to the other models [Fig. 5(C) and 5(D)]. When using modes m = 0, 1, and 2 (i.e., M = 2), the backscattering from the soft viscous-elastic model matches the modal solution and the finite element results (Fig. 5) over the whole frequency range (1–260 kHz). In Love's model, the peaks and nulls at higher frequencies are absent. Even though the first three modes suffice to include all peaks and nulls in the frequency range for the given bubble radius, Anderson's model is plotted using the first 20 modes to enhance accuracy of TS estimation as well as ensure that no peaks/nulls are excluded. To minimize the computational cost, the first three modes are included in the viscous-elastic model and shows a good agreement with Anderson's model over the whole frequency range.

Different parameters of the model such as acoustic, elastic, and damping properties as well as physical dimensions of scatterers were changed (Fig. 6) and compared to a base model (line 0) to investigate the effect on the backscattered energy over a frequency range of 1–260 kHz.

As can be expected, reducing the swimbladder radius (line 1) shifted the first resonance as well as the other peaks and valleys to higher frequencies. Increasing the shear modulus (line 2) or increasing the elastic layer (i.e., swimbladder wall) thickness (line 3) shifted the first resonance to higher frequencies while the other peaks and valleys did not shift. Increasing (here doubling) the fish flesh thickness (line 4) had a minor effect on the peak amplitude of the main resonance while slightly changing the amplitude of peaks and valleys at higher frequencies. Increasing (here tripling) the viscosity of flesh (line 5) decreased the main resonance amplitude and also had minor effects on the higher frequency peaks and valleys. Increasing the density of the gas inside the swimbladder (line 6), which corresponds to the fish being deeper, shifted all the peaks and valleys to higher frequencies. Furthermore, it widens the higher frequency peak and nulls and makes them more visible in the broadband backscattering measurements with their limited frequency resolution. The main resonance shifted more than the higher peaks and valleys. In addition to the gas density, increasing the flesh viscosity (line 7) will reduce the amplitude of resonance. For the increased density compared to the base model, reducing sound speed of the gas inside the swimbladder (line 8) by 10% (e.g., a carbon dioxide and oxygen mixture instead of pure oxygen), slightly shifts the main resonance to a lower frequency but with a more significant reduction in the higher frequency peaks and valleys. Bulk viscosity had negligible effect on the backscattering (not shown).

To select the model parameters systematically and to fit the model to data, the following protocol was used:

1. The density and sound velocity of oxygen inside the swimbladder, and in seawater were estimated from Eqs. (A13)–(A15) and Massel (2015), Appendix A, using the measured in situ temperature and salinity at the depth where a given target was present.

2. Shear modulus and swimbladder wall thickness were selected based on available or assumed values from the literature (Tables IV and V and Sec. II B 2 d).

3. $R4$ (radius of swimbladder) was manually selected by trial and error to match the peaks and overall magnitude of the model to the measured TS. Afterwards, to have the fish neutrally buoyant, $R2$ was calculated [Eq. (18)]. The value of $R2$ depended on the assumed flesh density $ρ2$⁠.

4. $μ2$ (shear viscosity, fish flesh) was selected to adjust the magnitude of the peaks from model to the measured TS.

For all 12 manually selected targets there was good agreement between measured and modeled target strengths over the measured frequency range (38 kHz narrowband and 50–260 kHz broadband) (Fig. 7). To model the TS, the following parameters are fixed (Table IV): swimbladder wall thickness (⁠$Δ=0.1$ mm), shear modulus (⁠$μ3=1.0$ MPa) and Lamé's first parameter constants $λ=2.4$ GPa. The sound speed and density of gas and seawater are functions of environmental parameters and depth and were estimated per Secs. II B 2 a and II B 2 b. Subsequently, swimbladder radii (⁠$R4$⁠), thickness of fish flesh (⁠$R2$⁠) and fish flesh shear viscosity (⁠$μ2$⁠) were adjusted to fit the measured TS and are given in Table VI.

The estimated swimbladder radii (⁠$R4$⁠) for the twelve targets ranged from 0.29 to 0.80 mm (Table VI).

For the modeling results (Fig. 7), shear modulus of the swimbladder wall was assumed to be 1 MPa, which corresponds to the shear modulus of medium rubber or latex (Table IV) and is within previously reported values for the shear modulus of a swimbladder wall (Sand and Hawkins, 1973; Feuillade and Nero, 1998; Fine et al., 2016). The thickness of the swimbladder wall (Δ) was set to 0.1 mm (Table VI) which is within the range found for mesopelagic fish species (Marshall, 1960) (Table V). To investigate the effect of shear modulus and swimbladder wall thickness on target strength frequency response, we applied different values for the shear modulus (⁠$μ3=0.2, 1, 2$ MPa) and wall thickness (Δ) (Fig. 8).

For each set of assumed shear rigidity and wall thickness of swimbladder, swimbladder radius, $R4$⁠, and the shear viscosity, $μ$⁠, were selected (shown on Fig. 8) to manually fit the model to the measured TS (see Sec. III C). It was observed that $R4$ (swimbladder radius) is the same in all of the curves regardless of shear rigidity or wall thickness of swimbladder. On the other hand, the shear viscosity (⁠$μ2$⁠) and wall thickness $(Δ)$ are connected and hard to disentangle. By increasing the swimbladder wall thickness (⁠$Δ)$⁠, a larger value of the fish flesh viscosity, $ζ2$⁠, is required. This is probably because the viscous layer thickness is decreased (as $R4$ remains constant) and the required damping is thereby provided by increasing the viscosity of the flesh. A larger shear viscosity will reduce the amplitude of the resonance peak due to the damping effect but will lead to a negligible increase in the resonance frequency (Feuillade and Nero, 1998), which is also evident in Fig. 8.

To demonstrate that the model was able to explain the measured target strength frequency response data, we chose twelve representative targets (Fig. 7). However, not all the measured targets showed such a good correspondence with the modeled frequency response. Three representative examples for which the viscous-elastic spherical model did not provide a good fit to the measured TS of individuals are outlined in Fig. 9. For targets 13 and 14, the model overall fits well to the measured TS until about 150 kHz, where the modeled TS becomes lower and higher for target 13 and 14, respectively. Also, the model fits reasonably well around the main resonance of the measured TS but the higher frequency peaks and valleys do not. For target 15, the model does not well capture the measured TS over a wide frequency range.

In the model, the morphology of the swimbladder was a sphere. To investigate the effect of swimbladder shape on the target strength frequency response, backscattering from a gas-filled prolate spheroid was estimated using the finite element method (see Sec. II C). The sphere in that model is replaced by a prolate spheroid and its backscattering is estimated. Target strengths for a prolate spheroid with aspect ratio (major divided by minor axis length) of 2 and an equivalent spherical radius (ESR) of 1 mm were estimated at six incident angles (Fig. 10) together with TS from a sphere of radius 1 mm for comparison. Backscattering from a sphere is independent of incident angle.

When comparing the TS curves from a prolate spheroid and spherical swimbladder (Fig. 10), it is evident that the main resonance of an elongated spheroid occurs at a slightly higher frequency than for a sphere [Fig. 10(A)]. Peaks and valleys at frequencies above 50 kHz occur at a lower frequency for a spheroid [Fig. 10(B)] than those from a sphere [Fig. 10(C)], which could explain the measured TS for target 14 and 15 (Fig. 9).

It is furthermore evident that the incident angle is important with regard to the overall backscatter intensity at frequencies greater than about 70 kHz, where an incident angle of 90° gives rise to an increase in TS of up to approximately 6 dB compared to an incident angle of 15° (Fig. 10), which could explain the deviation between measured and modeled TS for targets 13–15 (Fig. 9).

The two-layer viscous-elastic spherical swimbladder model that we used to describe the wideband in situ measured TS from individual mesopelagic organisms, has twelve tunable parameters (Table III). The sound speed and density of seawater and gas inside the swimbladder are derived from in situ measured environmental variables and are estimated by existing equations. The other parameters are confined to be within the reported values in the literature. The resulting swimbladder size estimates (Table VI) are within the range of previously reported values for mesopelagic fish (Table VII). However, fish are not the only type of organisms that potentially have gas inclusions (and thereby resonance) at mesopelagic depths, as siphonophores with pneumatophores (Barham, 1963) also share this characteristic. However, based on the trawling data and observations, the measured gas-bearing targets are unlikely to be siphonophores. A previous study by Knutsen et al. (2018) observed high densities of siphonophores in some Norwegian fjords and after a haul, the trawl meshes had red “slime” (“fouling”) attached, which was identified partly as remnants of physonect siphonophores (Knutsen et al., 2018). Similar “fouling” was never observed in the present study. We did not have dedicated sampling for siphonophores, but data from a Multinet and a macroplankton trawl, together with optical sensors, were used to evaluate their presence at the station. No larger physonect siphonophores, or remains thereof, were registered in the catches from a 1 m² Multinet deployed at the station. A stereo camera system (Deep Vision, Rosen and Holst, 2013) was attached to the aft of the macroplankton trawl and took pictures of organisms entering the cod-end. No pictures of physonect siphonophores were identified at the station. One single individual, which we would identify as a physonect siphonophore, was observed in the Deep Vision images at∼500 m depth at another station than the one used here (data not shown). A video plankton recorder (VPR) (Davis et al., 1992) attached to the front of the towed body MESSOR, continuously obtained images in undisturbed water ahead of MESSOR during operation. While the VPR imaged a small volume, the effective detection volume for animals utilizing extensive capture nets, such as physonect siphonophores, is much bigger than from a fish since the capture nets typically cover large volumes. In the VPR data, no physonect siphonophores were positively identified in the depth region 300–900 m at the station, and in a total of 142 597 images captured by the VPR in this depth range, the images contained two extended capture nets that could potentially belong to siphonophores, whether they are gas bearing or not. In summary, the underwater imagery data and the physical catches suggest that siphonophores are not likely candidates for most of the echoes, and therefore we deduce that the gas-bearing targets detected within the mesopelagic zone are dominated by gas-bearing fish.

Our results suggest that the modeling scheme used in this paper is appropriate for representing the broadband backscattering of gas-bearing mesopelagic organisms. Other models which are based on monopole scattering from a sphere (Love, 1978) or prolate spheroid (Ye, 1997), are widely used in studying backscattering of mesopelagic organisms (Kloser et al., 2002; Scoulding et al., 2015; Verma et al., 2017) near the main resonance frequency. Since these models do not include the higher modes of scattering, they fail to describe some of the features in the higher frequency regions, and therefore the target strength information within measured high frequencies remains unexploited. In addition, swimbladder wall thickness, its shear modulus and fish flesh viscosity are estimated. However, these parameters are somehow connected and further refining one of them would result in reduced uncertainties of the other ones.

Despite the strength of the model used in the present paper to explain observed characteristic features of mesopelagic fish backscattering, it is indisputably an approximation and has limitations. One of the reasons is that information on the acoustic and mechanical properties of live tissues are scarce, especially for mesopelagic organisms. Another reason is that the swimbladder shape and structure is more complex than our model. The model did not succeed in describing all the measured target strengths frequency responses, as demonstrated in Fig. 9, which could indicate that a spherical model is not a good representative of their swimbladder or some of the model parameters are incorrect for those specific targets. The swimbladder shape is one feature demonstrated to affect the TS frequency response (Fig. 10). Measurements of swimbladders in some mesopelagic fishes have shown that the shape is more prolate spheroid-like than spherical [e.g., Marshall (1960) and Kleckner and Gibbs (1972)]. Correct backscattering estimation over a wide frequency range for a bubble with a prolate spheroid shape is more complicated than for a sphere especially for large aspect ratios (Prario et al., 2015). Analytic backscattering modeling of a gas-filled prolate spheroid including viscous and elastic layers would be highly computationally demanding. We are aware that our model is simplified in the sense of swimbladder shape, but have instead focused on including more physical parameters such as the damping effect of the fish flesh and swimbladder wall. To better understand the limitations of spherical swimbladders and to avoid over-interpretation, we modeled backscattering from a prolate spheroid with aspect ratio of 2 for different incident angles over a wide frequency range (Fig. 10). It was observed, and is well-known that the incident angle is important for the frequencies well above the resonance (Scoulding et al., 2015). Comparing the backscattering of spherical bubbles, spheroidal bubbles, and viscous-elastic spherical gas bubbles of the same volume over a wide frequency range (Figs. 6 and 10) would help to understand the significance of each parameter in different frequency regions.

Previous studies applying a prolate spheroid resonance model to mesopelagic gas-bearing organisms have found that the modeled TS was sensitive to swimbladder volume and aspect ratio, tilt angle, and viscosity of the fish flesh (Scoulding et al., 2015; Proud et al., 2019), which in turn can result in uncertainties when estimating biomass based on backscattering energy. The model in the present study can (ideally) be applied in order to obtain information on physical characteristics of the swimbladder, such as shape and radius and thereby decrease uncertainties in biomass estimations of mesopelagic fish. In most fish species, the swimbladder size is positively correlated with size of the fish, which enables estimates of fish size based on their TS at a given acoustic frequency [e.g., Nakken and Olsen (1977)]. However, for some mesopelagic fish species, the volume of the swimbladder has been found to decrease with size (Marshall, 1960; Butler and Pearcy, 1972; Davison, 2011; Scoulding et al., 2015), thereby complicating size estimation based on TS measurements. This emphasizes the importance of looking further into the intraspecific and interspecific differences of mesopelagic fish physical characteristics. This was also noted by Scoulding et al. (2015) who found different swimbladder anatomies (fish length dependence and swimbladder presence/absence) for two mesopelagic fish species belonging to two different families: Myctophidae and Sternoptychidae.

Different gas content/composition inside the swimbladder will change the acoustic properties (density and sound speed). In the viscous-elastic model it was assumed that the swimbladder was filled with oxygen. The density and sound speed inside the swimbladder was then estimated per Sec. II B 2 a. The gases present in a swimbladder are the same as gases of air dissolved in water [Fänge (1966), and references therein], and the gases have different sound speed properties where some have higher and other lower sound speed in comparison with oxygen. Changes in the sound speed of gas inside the swimbladder will have a larger effect on the higher frequency peaks and valleys than on the main resonance. For example, by reducing the sound speed of a gas while keeping the density unchanged, the higher peaks and valleys shift to lower frequencies (as observed by comparing curves 6 and 8 in Fig. 6). Accordingly, the model has a better fit to the measured TS of target 13 (Fig. 9) if the swimbladder is filled with a gas with a higher sound speed than oxygen. Different gases pass through biological membranes at different rates (Krogh, 1919) and thus, the composition of gas in a swimbladder will differ with depth and between individuals and species (Fänge, 1966), which could explain the deviation of the model from the measured data in some cases. Also, some mesopelagic fish species have fat-invested swimbladders (Marshall, 1960; Butler and Pearcy, 1972), which in turn could have an effect on the TS as the density and sound speed inside the swimbladder would change.

For some of the in situ measured and modeled target strength frequency responses in the present study, some features (TS peaks and nulls) were observed at higher frequencies (Fig. 7). These features could possibly provide information about morphological characteristics such as length and width of the fish and shape of swimbladder [e.g., Reeder et al. (2004) and Kubilius et al. (2020)]. For example, peaks and nulls could be consistent with constructive and destructive interference between different parts of the body as observed in the swimbladder-bearing fish species alewife (Alosa pseudoharengus) (Reeder et al., 2004). If the target strength frequency response could provide both the flesh weight and size of the swimbladder, this has the potential to enable identification and separation of different species and size groups in mixed scattering layers. Thus, by tuning the different parameters in the model presented here, it might be possible to estimate some distinct physical and biological characteristics such as the ones mentioned above.

Calibration of echosounders are conducted mainly to calculate the on-axis gain and the equivalent beam angle (Simmonds and MacLennan, 2005). The in situ measurements were collected with echosounders calibrated at the surface and not at the measurement depth. Temperature and pressure influences the performance parameters of echosounders (Kloser, 1996; O'Driscoll et al., 2013; Haris et al., 2018). However, older-technology echosounders with air-back design transducers are more affected than the newer ones with composite transducers (Haris et al., 2018), which are similar to the ones (70, 120, 200 kHz transducers) used in the present study. Calibration from surface to 900 m depth affects the on-axis gain for the composite transducers by less than 1 dB (Haris et al., 2018). The equivalent beam angle at depth is less than the nominal values given by the manufacturer (Haris et al., 2018). However, the latter has only minor effects on the results in the present paper, since we only used targets from near on-axis (within 3°).

In this paper, the procedure and applicability of the method were demonstrated by manual single target detection and fitting of a viscous-elastic model to their measured TS. To apply the demonstrated procedure to large datasets, effective automated single target detection algorithms for broadband acoustic data and an optimization algorithm to fit the model to the measured TS of detected single targets, are required.

Net data and optical observations revealed very low abundances of siphonophores at the station, which lead us to presume the acoustic targets were dominated by gas-bearing fish. We used a viscous-elastic spherical swimbladder model which includes higher modes of backscattering, and hence could explain the observed features of the in situ measured backscattering from individual mesopelagic gas-bearing fish. Since the model is physics-based, the modeling parameters can be related to the physical properties of fish. By combining broadband acoustic target strength measurements with the mathematical/physical backscattering model, it was possible to obtain physical characteristics of the targets such as swimbladder size. The model can further be applied to obtain additional physical information about the targets by investigating the other tunable variables resulting in, for example, information on swimbladder shape, which is one parameter leading to uncertainties in biomass estimations of fish based on backscatter measurements. By using broadband acoustic backscatter, some detailed features were observed in the frequency responses, which were also successfully captured by the model. In addition to this, applying the model to measured backscattering data could in turn be used to estimate the flesh weight of the fish, assuming neutral buoyancy. The proof of concept demonstrated in this paper should in the future be automated and applied to large datasets. This would require more effective single target detection algorithms for broadband acoustic data and an optimization algorithm to automatically fit the model to the measured TS of detected single targets.

We thank the officers and crew onboard R/V Kronprins Haakon. Tonje Nesse Forland is thanked for helping with the FEM backscattering modeling. Benjamin Marum is thanked for helping with data collection. Funding for this work was provided by the Institute of Marine Research (project No. 15093), HARMES project, Research Council of Norway (project No. 280546), and MEESO, EU H2020 research and innovation programme (Grant Agreement No. 817669).

The velocities and stresses based on the velocity potentials in the spherical coordinates (Karlsen and Bruus, 2015),

Stresses for the fluid medium,

Stresses for the solid medium,

The Legendre differential equation in spherical coordinate solved by Legendre polynomials $Pm(cos θ)$ (Arfken and Weber, 2005),

Some of the boundary conditions. The first boundary condition (⁠$vr1r=R2=vr2r=R2$⁠) will be

The fifth boundary condition (⁠$σrr2r=R3=σrr3r=R3$⁠) will be

The seventh boundary condition (⁠$vθ2r=R3=vθ3r=R3$⁠),

and the tenth boundary condition (⁠$σθr3r=R4=0$⁠)

Thermodynamic equations for density and sound speed in fluids. The ideal gas equation of state is (Boyle's law)

and Van der Waals equation of state is