A thorough understanding of the scattering characteristics of marine organisms is a prerequisite for robust quantitative fisheries acoustic data processing or interpretation. Target strength models, such as the distorted wave Born approximation (DWBA) can be used to improve the understanding of field recordings of weakly scattering targets. With acoustic methods now being used by a wide audience, allowing access to such models becomes a necessity. To ease access to the DWBA model, an r package (zooscatr) which includes a web application and the ability to parameterise the model either through the web application, text files, or pure scripting has been developed and is now freely available on Github.

## 1. Introduction

The use of acoustic technologies as a standard tool for the monitoring and characterisation of marine organisms requires a detailed understanding of the sound scattering signatures of organisms within marine ecosystems. With the emergence of broadband acoustics, acoustic scattering models can help us to better understand the received acoustic signal at an improved resolution, when compared to narrowband acoustics. Among the scattering marine organisms, weakly scattering (fluid-like) zooplankton or crustacean species, such as euphausiids, salps, and copepods, make up a significant portion of the global marine biomass (Strömberg *et al.*, 2009). To accurately estimate their biomass acoustically, it is essential to provide accurate acoustic scattering models that allow for fluctuations in acoustic scattering as a function of their diverse morphological, physiological, and behavioural aspects, and acoustic frequency (Amakasu *et al.*, 2017; Bassett *et al.*, 2016; Stanton and Chu, 2000). The availability of fast, standardised, and flexible sound scattering models, accessible to a broader audience hence becomes desirable.

Jech *et al.* (2015) discussed the usefulness of the distorted wave Born approximation (DWBA) and other scattering models in fisheries acoustics for a number of targeted organisms and compared model outcomes to analytical solutions. The DWBA has widely been used to model the scattering properties of weak scattering targets, such as zooplankton and small crustaceans, e.g., krill or decapod shrimp (Chu *et al.*, 1993; Chu and Ye, 1999; Stanton *et al.*, 1996, 1998; Stanton *et al.*, 2000; Stanton and Chu, 2000). Based on previously developed matlab scripts (Chu, 1999), we developed an r (R Core Team, 2018) package called zooscatr, that allows for parameterisation through .dat or .csv plain text files, a web application (Fig. 1) or r scripting. zooscatr can deal with arbitrary cylindrically symmetric elongated shapes, with an arbitrarily varying radius along the cylindrical axis. Material properties can be variable along-axis and different lengths as well as target orientations can be defined. zooscatr is open-source and open-access. zooscatr (v 0.3.0) is provided in the supplementary material^{1} with any future developments available through Github (Gastauer, 2019).

## 2. The DWBA model

A detailed description of the applied phase-compensated DWBA (PC-DWBA) model and the application to scattering properties of zooplankton can be found in Stanton *et al.* (1998) and Chu and Ye (1999). This variation of the DWBA model was chosen because parameterization of the model is flexible (not restrictive in terms of scattering geometry, target shape, heterogeneous material properties, and acoustic frequency), computational speed and the resulting applicability to a wide range of weakly scattering organisms. Further, the model has been validated within the context of previous studies (e.g., Amakasu *et al.*, 2017; Chu *et al.*, 1993; Chu *et al.*, 2000; Chu *et al.*, 2003; Stanton *et al.*, 2000; Stanton and Chu, 2000; Wiebe *et al.*, 2009). Limitations of the DWBA approach have been discussed previously (Chu and Ye, 1999; Jech *et al.*, 2015). The DWBA model is only valid for weak scatterers, with material properties (sound speed and density contrast) within 5% variation of the surrounding medium and shear waves are assumed to be absent (Stanton *et al.*, 2000). Different variations of the DWBA have been proposed (i.e., Calise and Skaret, 2011; Demer and Conti, 2003; Jones *et al.*, 2009), that provide comparable results with minor variations mainly caused by differences in how the model deals with phase and shape parametrisations. The DWBA for a body of an arbitrary shape and size for a given frequency can be written as (Chu *et al.*, 1993)

where *f _{bs}* is the scattering amplitude;

*k*is the wavenumber (defined as

*2π/λ*, where

*λ = f/c*is the acoustic wavelength [m], with frequency

*f*[Hz] and

*c*the sound speed [m/s] in the surrounding fluid); the subscript

*s*refers to the surrounding fluid, the subscript

*b*refers to the scattering body; $rv\u2192$ is the position vector along the body axis;

*v*is the integration volume of the body; $\gamma \kappa $ and γ

_{ρ}are related to the body material properties [Chu

*et al.*, 1993; Stanton

*et al.*, 1998; also see Eq. (4)]. The relationship between the backscattering amplitude, the normalised backscattering amplitude (

*f*[ratio]), the scattering pressure (

_{norm}*p*), the differential backscattering cross-section (

_{scat}*σ*[m

_{bs}^{2}]), target strength (

*TS*[dB re m

^{2}]) and the reduced target strength (

*RTS*[ratio]) are summarised in Table 1 in the supplementary material.

^{1}

Within zooscatr, the DWBA model can be run through the *bscat*() function or the web application *DWBAapp*() (Fig. 1), which relies on *bscat*(). The model takes four categories of input parameters: shape (Table 2 in the supplementary material^{1}), orientation (Table 3 in the supplementary material^{1}), material properties (Table 4 in the supplementary material^{1}), and simulation parameters (Table 5 in the supplementary material^{1}). Templates for the input parameter files are included in the package and detailed instructions on how to run the model can be found in the package vignettes.

The modelling function runs through three phases: (1) the shape generation, where the coordinates are defined; (2) the modelling phase, where all other parameters are read into the model and the scattering model is computed; and (3) the conversion phase, where the backscattering amplitude is translated into the selected output format (Table 1 in the supplementary material^{1}).

## 3. Shape generation

All shape inputs of the scattering objects are required to have circular cross sections along the length axis [shapes are formed using the *buildpos* (build positions) function]. Shape information contains *n _{i}* points along the body axis (defined within the simulation parameters), normalized

*x*and

*z*coordinates, tapering information (i.e.,

*y*coordinates), the positional vector $rpos\u2192$, a distance vector

*dr*, the angle

_{1}*γ*[rad] which equals the angle of a hypothetical line passing through the coordinates and the origin as well as the angular position of the points

_{t}*θ*(Fig. 2). zooscatr can be used to either model uniformly bent cylinders or arbitrary shapes made up of a series of adjacent, potentially tapered, cylinders.

_{tilt}A uniformly bent cylinder with regular tapering is the default shape computed, if a shape input is not given. Construction of a uniformly bent cylinder with regular tapering requires the ratio of the radius of curvature of the longitudinal axis (*ρ _{c}*) over

*L*(

*ρ*), the tapering order (

_{c}/L*n*) [Eq. (3)] and the ratio of

*L*[mm] over the cylindrical radius, i.e., half width of the cylindrical cross-section (

*a*[mm]) (

*L/a*) (Table 2 in the supplementary material

^{1}). Arbitrary shapes are specified by

*x*,

*z*position (Fig. 2) and tapering information and optionally

*x*and

_{upper}= x + tapering*x*.

_{lower}= x − taperingWhen an arbitrary shape is used, the coordinates for the *n _{i}* points are normalised with the range of

*x*and

*z*coordinates minus the range average. The length

*L*of the shape is computed as the cumulative sum of the two-dimensional distance of the coordinates

*x*and

*z*,

The final length normalised coordinates *x* and *z* are defined as $\u20092x/L$ and $\u20092z/L$, respectively, to fit within a range of 0 and 1. The tapering at each integration point is computed through linear interpolation of the given tapering points along the *z*-axis. Tapering is controlled by the tapering order (*n*) where the radius for a given point is defined as

Note that *n* = 2 corresponds to a prolate spheroid. A smoother function can be applied to both, the body axis (shape complexity) and tapering. The smoother function is a simple symmetrical moving average, where the smoother factor *n _{s}* will be set to the number of input values (

*n*) if it is larger than

_{l}*n*.

_{l}Independent of the shape, the length of the positional vector $rpos\u2192$ is computed as $x2+z2\u2009$. The angle between the origin and the coordinates in radians is calculated as the two-argument arctangent of the coordinates [$\gamma t=arctan2(z,x)$] (Fig. 2). The two-argument arctangent of the coordinates *z*, *x* represents the real and imaginary components of the complex number *z + ix*, allowing the user to define a unique phase in four quadrants (Fig. 2). The two components of the directional vector $dr\u2192$, *dx* and *dz*, are defined as $dx=xn+1\u2212xn$ and $dz=zn+1\u2212zn$.

## 4. Modelling phase

Besides orientation and length, the material properties, i.e., the sound speed contrast (*h*) and the density contrast (*g*) are important model input variables. These contrasts are explicitly contained within the general DWBA equation under the form of the term *C _{b}*,

$\gamma \kappa \u2009$ is related to the compressibility $\kappa =1/\rho c2$ and $\gamma \rho $ is defined through *g.* Both contrasts are the ratio of the material properties inside the body and the surrounding fluid, i.e., density (*g = ρ _{b}/ρ_{s}*) and sound speed (

*h = c*), respectively. Ideally

_{b}/c_{s}*g*and

*h*should be estimated as accurately as possible for the given target at ambient conditions (Chu and Wiebe, 2005; Matsukura

*et al.*, 2009a). Further, for euphausiids, a life history dependency as well as seasonality of

*g*and

*h*has been reported (Chu and Wiebe, 2005; Køgeler

*et al.*, 1987; Matsukura

*et al.*, 2009b; Mikami

*et al.*, 2000).

*g*can be measured by using volume displacement density bottles (Foote, 1990; Greenlaw, 1979; Mukai

*et al.*, 2004), a density increment method (Becker and Warren, 2014, 2015), or through dual density measures (Chu and Wiebe, 2005). The establishment of

*h*can be determined either through a velocimeter probe (Chu and Wiebe, 2005; Foote, 1990; Matsukura

*et al.*, 2009a) or through computational tomography (CT) scans of the targets (Au and Fay, 2012; Gastauer

*et al.*, 2017). If

*g*and

*h*cannot be determined directly, literature values should be consulted and model estimates should be surrounded by an appropriate uncertainty estimate.

If a homogenous body along-axis is assumed, density (*g*) and sound speed (*h*) contrasts are directly read from the parameters. For the case of an inhomogeneous body along-axis, the input *g* and *h* values are treated as means where the number of segments and the standard deviation of *g* and *h* must be defined. If the number of segments is greater than 8, the correlation length (% of body length) must be provided as well. The correlation length defines a larger variability around the midsection, with the variability decreasing towards the ends. This correlation length is the standard deviation of a Gaussian function centred at the midsection of the body.

Solving the DWBA equation with the provided parameters for an elongated fluid object with a circular cross section leads to the following equation for each orientation angle *j*:

With *ka* the wavenumber (*k*)∗ cylindrical radius (*a*), *t* the tapering coefficient, and *J*_{1} the Bessel function of the first kind of order 1. Prior to transforming the model outcome into the selected output variable format, if selected, the length and angular averages will be computed. If no average should be computed, only the mean value will be considered. For angular averages a uniform or a Gaussian PDF can be calculated

With *m* = the number of orientation points within a user defined angle range, and $\sigma \theta =sd(\theta tilt)$ the standard deviation of *θ _{tilt}*. The resulting scattering amplitude

*f*will then be computed as the square-root of the conjugate product of

_{2}*f*multiplied by the PDF of the angular distribution as $f2=f\u2009f*\u2009PDF$ and the differential backscattering cross-section becomes $\sigma bs=f22$.

For length averages either a uniform or a Gaussian PDF can be computed, where only the mean length will be considered if no length average should be computed,

where $m=6[sd(L)/dL]$ with *dL* the length increasing step, and $\sigma L=sd(L)$ the standard deviation of *L* After the length PDF has been defined, the updated *σ _{bs}* becomes the product of

*L*

^{2}, the length

*σ*over

_{bs}*ka*. The frequency-dependent differential backscattering cross section averaged over angle of orientation and length, $\u27e8\sigma bs\u27e9$ can then be calculated as $\Sigma \sigma bsi/m$. Following Table 1 in the supplementary material,

^{1}the target strength or reduced target strength can then be deduced.

## 5. Model validation

Following methods described in Jech *et al.* (2015) the analytical solution for a weakly scattering fluid sphere based on Anderson (1950) and for a hypothetical fluid sphere with the material properties of a generic krill (Calise *et al.*, 2011), ensonified by a plane wave, were used to benchmark the model. A summary of the sphere parameters can be found in Table 6 in the supplementary material.^{1} A detailed description of the benchmarking process is included in zooscatr as a vignette.

The overall mean error was computed following Jech *et al.* (2015),

The overall mean error for the zooscatr implementation was 0.31 dB, which is comparable to results reported by Jech *et al.* (2015) for the DWBA model (0.38 dB) [Fig. 3(a)]. While the model performed very well in general, main variations to the analytical solution occurred at the location of the nulls, where the zooscatr model generally delivered lower TS values (Fig. 3). Similarly, when comparing the model output to the analytical solution for a fluid sphere with material properties of krill [Fig. 3(c)], the model generally performed well, detecting the troughs and peaks at the correct locations [Fig. 3(d)]. The overall average error was 1.93 dB, but if only TS values > −80 dB were considered, the error was reduced to 0.76 dB. Lower values, detected within the null locations, are considered to be of low importance, as these are likely to be linked to precision and the resolution at which the model is run rather than an error in the model itself.

## 6. Conclusions and outlook

Comparing modelling results to analytical solutions has shown that the DWBA model holds validity for weakly scattering targets. The next logical step would be the validation of the model with *in situ* or *ex situ* information. With the availability of more and more small and flexible acoustic sampling platforms with broadband capabilities, in combination with optical recordings, the availability of fast and flexible models could help to improve our understanding of the backscattering biomass and behaviour of weak-scattering layers found throughout the world's oceans. At the moment of writing, the inversion function is a simple lookup function, minimising a cost function defined as a sum of squares error function. This is a very basic function, which should be treated with care and seen as a first step towards more sophisticated inversion methods. A benefit of the provided implementation of the DWBA model within r is its capability of running the model at relatively high computation speeds, allowing the simulation of millions of targets within minutes (on a standard notebook with six cores). Inversion accuracy could be enhanced through the application of machine learning (MLA) and deep learning algorithms (DLA). The results of large numbers (millions) of DWBA simulations will provide information gain on the variability of the acoustic scattering properties of weakly scattering targets, in respect to their shape, orientation, internal structure and the selected frequencies and could form the basis for robust MLA and DLA training datasets.

Compared to previously developed tools [such as the online DWBA calculator developed by the NOAA Fisheries Southwest Fisheries Science Center (https://swfscdata.nmfs.noaa.gov/AST/SDWBA/sdwba.html)], major benefits of the present package are that it is open-source and open-access, providing the DWBA modelling capabilities to a wider audience. In order to make the package more accessible to the users, several vignettes have been included in the package, which can be read as tutorials, describing, among other things, how to use the integrated web application or how to use it from within the r command, how to define a parameter file, and how to run the model in parallel mode for the fast simulation of numerous targets.

Future developments of the package may include more sophisticated inversion methods such as an acoustic scattering model-based nonlinear inversion algorithm (Chu *et al.*, 2016) or MLA/DLA algorithms, the inclusion of other scattering models, porting the model processing from the CPU to the GPU would increase the speed substantially.

^{1}

See supplementary material at https://doi.org/10.1121/1.5085655 E-JASMAN-145-504901 for the zooscatr package binary (v.0.3.0). Support and package updates will be provided through Github (Gastauer, 2019). The package is distributed under the GNU Public License as published by the Free Software Foundation, either version 3 of the license, or (at the reader's option) any later version. A definition of the scattering metrics used within ZooScatR and their relationship to the scattering amplitude is given in Table 1. A description of the model shape input (Table 2), orientation (Table 3), material properties (Table 4), and simulation (Table 5) parameters are provided. The used model input parameters for a weakly scattering fluid sphere can be found in supplementary materials (Table 6).