Modelling the three‐dimensional space use of aquatic animals combining topography and Eulerian telemetry data

Passive acoustic telemetry provides the opportunity to monitor and contextualize the movements of diverse aquatic animals. Despite depth being an essential dimension along which many processes are organized, the Eulerian structure of the acoustic telemetry data (movements perceived from fixed locations) and the consequences of sound propagation in water hinders the incorporation of the vertical dimension into animal’s space use analyses. Here, we propose a new data‐driven quantitative method to estimate 3D space use from telemetry networks. The methodology is based on simulating large numbers of stochastic synthetic paths, accommodating the detection probability around receivers and the depth information from transmitters and integrating the local topography. The methodological protocol is explained in detail and tested with acoustic telemetry data from the common dentex Dentex dentex in a Mediterranean marine protected area. We present 3D space use estimations for the tagged individuals and compare them with other 3D and 2D estimations derived with existing probabilistic methods. 3D space use estimations that incorporate topography provided a more comprehensive view of the movement ecology of tracked individuals, with relevant pieces being missed by 2D representations. Our method generated realistic representations of the actual spatial co‐occurrence of individuals, including the spatio‐temporal identification of relevant aggregation areas.

Acoustic telemetry is the most widely used technique to monitor the movement of aquatic animals (Hussey et al., 2015). In passive acoustic telemetry, animals tagged with acoustic transmitters are monitored by a network of receivers placed at fixed locations in the study area (Heupel, Semmens, & Hobday, 2006), and thus, movement is perceived in an Eulerian mode as presences, absences or transitions between discrete locations. The uncertainty regarding the position of tagged animals is usually large (hundreds of meters), as it depends on the distance at which acoustic signals are detectable (i.e. the detection range).
Accurate depth data are often provided by the acoustic transmitters attached to animals (Veilleux et al., 2016), but this information is rarely taken into account when estimating the space use of fish. In recent years, several probabilistic approaches have been proposed to quantify volumetric space use by incorporating the vertical dimension into location-based methods, such as kernel density estimations (Simpfendorfer et al., 2012;Worton, 1989), or into movement-based methods, such as the Brownian bridge movement model (Horne, Garton, Krone, & Lewis, 2007;Kranstauber, Kays, Lapoint, Wikelski, & Safi, 2012;Tracey et al., 2014). However, these methods require a higher precision in the locations of animals, which is only possible when acoustic receivers are placed on a relatively close grid, allowing the calculation of centres of activity (Simpfendorfer, Heupel, & Hueter, 2002) or the application of more sophisticated triangulation methods (e.g. VEMCO Positioning System;Steel, Coates, Hearn, & Klimley, 2014). Typical coast-lined configurations of acoustic arrays are problematic given that non-mesh architectures prevent highly resolved spatial locations and complicate the aforementioned volumetric space use estimations.
Another challenge of movement analyses in two and three dimensions is the incorporation of topographical characteristics affecting space use estimations. Prominent types of landforms, such as mountains or rivers in terrestrial ecosystems, or emerged landmasses in aquatic ones can pose insurmountable barriers to the movement of animals, a fact that should be incorporated in the estimation of the space use, yet rarely is included in the analysis (e.g. Benhamou & Cornélis, 2010;Wilson, Hanks, & Johnson, 2017).
Here, we propose a new computational method to estimate 3D space use probability distributions from passive acoustic telemetry data. It is an empirically driven method, which integrates the location uncertainty of coarse acoustic telemetry data by simulating large numbers of stochastic 3D trajectories (i.e. synthetic paths) according to the observed detection probabilities around receivers and the movements of individuals within the array.
Our method assumes simple movement characteristics to generate synthetic paths, which follow the shortest pathway between two points and take into account the topography of the study site to avoid movements through hard boundaries (i.e. emerging landmasses and shallow areas). This contrasts with other 2D and 3D space use approximations grounded on the Brownian motion model, which assume a diffusive movement model defined by a single parameter, the Brownian motion variance (Horne et al., 2007;Kranstauber et al., 2012;Tracey et al., 2014). Our method, instead, incorporates the variability observed in the original data and constructs consistent space use estimations by integrating large numbers of iterations of synthetic paths. This approximation accounts for movement barriers in a natural way, but implies the use of contextual data (physical characterizations of the study site) beyond the one provided by the acoustic array.

| OVERVIE W OF THE WORKFLOW
The method to estimate 3D space use probabilities has been developed in the r programming environment (R Core Team, 2017) and implemented in the fishtrack3d package (https ://github.com/aspil laga/fisht rack3d). The workflow is divided into three steps: (a) empirical characterization of detection probabilities around receivers, (b) generation of synthetic paths and (c) assembly of 3D volumes.
The main features of each of the steps are sketched below, but we provide a detailed explanation in a vignette included as supporting information (Appendix S1) and within the fishtrack3d package.

| Empirical characterization of detection probabilities
Two main factors modulate the probability of being detected by an acoustic receiver: the distance between the transmitter and the receiver (in relation to the acoustic range) and the presence of physical impediments for the transmission of acoustic signals. The acoustic range is tested by placing several receivers and one or more transmitters at increasing distances in between to determine the ratio of emitted signals that are detected at each distance (Kessel et al., 2014). We model the average decay in detection probability with distance by fitting a logistic regression model to an empirical range test carried out in the study area. In addition, the presence of emerged landmasses or heterogeneous rocky bottoms generate areas of acoustic shadows within the potential range of the receiver array, from which tagged animals cannot be detected. We use a viewshed analysis to identify the acoustic shadows generated by large emerged landmasses, to then incorporate their effect during the simulation of synthetic paths. Viewshed is a commonly used terrain analysis technique that recognizes the geographical area that is visible from a specific location (Haverkort, Toma, & Zhuang, 2009), which is also applicable to sound propagation to detect, over a large scale, the areas from which acoustic transmitters could be detected by each receiver.
The acoustic range and the distribution of the acoustic shadows present spatial and temporal variations due to physical and biological factors, such as sound diffraction, hydrographic conditions (e.g. currents, swell), habitat characteristics and environmental noise (Heupel et al., 2006), which are omitted in this first version of our method due to the difficulty of quantifying them. However, one of the advantages of our method is its high flexibility to accommodate additional data and further modifications in future, such as the inclusion of habitat-or receiver-specific acoustic range models that account for spatial heterogeneity (e.g. using mixed effects models, Selby et al., 2016), or information from fixed control tags to account for temporal variation. Moreover, high-resolution cartographies and voxel-based 3D viewshed algorithms might be used to improve the characterization of acoustic detection landscapes, at the expense of increasing the computational cost.

| Generation of synthetic paths
In the next step, we simulate a large number (e.g. 100) of stochastic trajectories (i.e. synthetic paths) from the telemetry data. Each synthetic path requires a computationally intense procedure, and its duration depends on the number of detections to be processed. Therefore, a prior thinning is highly recommended for large datasets with high temporal resolution (e.g. detection intervals of < 15 min) and long tracking periods. In order to reduce the amount of data while maintaining its original variability, we propose a random thinning method that pools telemetry data into specific time intervals (e.g. 30 min, 1 hr), generating a unique combination of data to start each simulation of a synthetic path (see Supporting Information Appendix S1 for more details). By adjusting the length of the intervals, this method can be also used to combine data from transmitters with different emission rates.
We construct the synthetic paths by sampling a pair of geographical coordinates for each detection associated with a receiver, considering the detection probability, the depth values from transmitters and the local topography (bathymetry). The coordinate sampling is based on a probability raster, obtained by applying the acoustic range model to the distances between a given receiver and the raster cells at the depth of the observation. In addition, we consider the distance to the previous location by also excluding the cells that are beyond a maximum distance that the animal might have reached in the time elapsed from the previous observation, assuming a fixed maximum speed (i.e. 1 m/s).
This allows us to avoid unlikely movements between distant locations in short periods of time. After the sampling process, consecutive pairs of 3D coordinates are joined using a least-cost path approach, which finds the shortest route between locations avoiding to cross unsuitable areas such as emerged landmasses or shallow cells (see Supporting Information Appendix S1 for more details). As a final result, we obtain a unique three-dimensional synthetic path from each iteration.

| Assemblage of utilization distribution volumes
Three-dimensional utilization distributions (3D-UDs) are finally estimated from the average spatial occurrence of all the synthetic paths. We present the study area as a three-dimensional grid of voxels with fixed horizontal and vertical resolution (i.e. 10·10·1 m 3 in our example), to then calculate the average percentage of time that synthetic paths spend in each voxel. These average values are finally processed by applying a three-dimensional kernel estimation, using the "ks" package for r (Duong, 2007), to generate smooth 3D-UDs that are easier to manipulate and plot. 3D-UDs can be used to calculate the volume sizes and spatial overlaps for different probability contours, or be plotted using 3D visualization engines such as the "rgl" (Adler et al., 2018) and "plotly" (Sievert, 2018) packages for r or external software such as paraview. (Ahrens, Geveci, & Law, 2005).

| TE S TING THE ME THOD
The method was first tested with a simulated acoustic telemetry dataset, corresponding to a set of individuals moving through a given receiver array following some simple pre-established parameters. We compared the outcome of the method, with and without including topography, with the results that would be expected from the pre-established parameters. The procedure to simulate the data and the results from this test are provided in supporting information (Supporting Information Appendix S2).
In addition, we applied and tested the performance of the method with real acoustic telemetry data from the common dentex as a surrogate of path lengths and increased the number of synthetic paths from 1 to 100. Finally, we quantified the spatial co-occurrence of the common dentex during the spawning period (taking place between May and June), when individuals punctually migrate to deep rocky bottoms to breed (Aspillaga, 2017). In particular, we computed the overlap between 3D-UD week-long estimations for some tracked individuals outside and within the spawning period, separately for day and night periods as defined by local sunset and sunrise times provided by the NOAA Solar Calculator (www.esrl.noaa.gov/gmd/grad/ solca lc/).

| Optimizing the number of synthetic paths
Generating each synthetic path proved a computationally intense and time-consuming step (45 ± 8 s for a 100-location path on a computer with 4GB RAM); hence, the method might be optimized by adjusting the number of synthetic paths that are required to reach a stable 3D-UD solution. Short tracking periods consisting of fewer locations require large numbers of synthetic paths in order to reduce the variability of the resulting space use estimations ( Figure 1). Long tracking periods, by contrast, contain many more data points that already capture most of the space use variability; hence, robust space use volumes are obtained with fewer synthetic paths. In our case, at a temporal resolution of 30 min, one-day-long tracks (41 ± 8 fixes per track) required almost 70 synthetic paths to stabilize the size of the space use volume estimations. In contrast, 30-or more-day long tracks (more than 1,382 ± 60 fixes per track) required only half of the sampling effort, that is about 35 synthetic paths (Figure 1).

| 3D versus 2D space use estimations
The 3D space use probability estimations for common dentex in- When considering the horizontal plane only, the results obtained from both the collapsed 3D-UDs and 2D-UDs calculated with the dBBMM were highly consistent (Figures 3 and 4b). This robustness is remarkable given that the collapsed 3D-UDs make minimal assumptions on the underlying nature of the movement process, while the 2D-UDs were calculated assuming a conditional Brownian movement between locations (Horne et al., 2007;Kranstauber et al., 2012).
The most relevant differences between the resulting 3D space use estimates were related to the topography, which allowed the effective exclusion of unfrequented shallow areas and unsuitable emerged landmasses and the redistribution of the UD in those areas that were in fact suitable for individuals (Figures 3 and 4c, and App #2 in Supporting Information Appendix S3). By handling and incorporating geographical barriers into space use characterizations, our data-driven approach clearly led to more accurate and realistic 3D-UD and 2D-UD estimates (Supporting Information Appendix S2).

| Spatial overlap during the spawning season
A 3D approach provides a unique view of the space use patterns of fishes, in particular, when studying the spatial overlap between individuals. The overlap indices of any given pair of 3D-UDs were lower on average than the overlap values of the 2D-UDs, but without any apparent correlation between them (Figure 4a). Correspondingly, the aggregating behaviour of common dentex during the spawning season (Aspillaga, 2017) is strongly evidenced when analysing 3D overlaps. At night, the 3D overlap between pairs of individuals was significantly higher in the spawning season compared to the non-spawning season (Kruskal-Wallis rank-sum test, H = 34.6, df = 1, p-value < 0.05, overlaps for the 95% contour level; Figure 5). These differences were not significant if 2D-UD overlap data were used (Kruskal-Wallis rank-sum test, H = 2.8, F I G U R E 1 Relative change of 3D-UD volumes for 50% (a) and 95% (b) contour levels when increasing the number of synthetic paths and their length (number of simulated successive days) F I G U R E 2 3D-UD volumes (50% contour level) of four common dentex individuals (different colours and codes) showing a depth segregation. Vertical axis is 10 times magnified compared to the horizontal axes. Violin plots on the right show the distribution of depth values of the original detections. The 3D-UD volumes of all the common dentex individuals can be interactively explored in the shiny app #1 in Appendix S3 F I G U R E 3 Bi-dimensional representations of utilization distributions of two common dentex individuals (rows) calculated with different methods: the dynamic Brownian bridge movement model (dBBMM, first column), 3D kernel density estimation on centres of activity (CoA, second column) and the new proposed method without and with incorporating topographic information (third and fourth columns, respectively). 3D-UDs were collapsed to the horizontal plane by vertically summing the UDs of all the depth layers. Polygons indicate, from the inside out, the 50% and 95% contour levels This example illustrates why topography-wise 3D-UD estimations are a valuable tool to detect and visualize susceptible areas, providing precise information to focus on conservation efforts.

| CON CLUS IONS
We presented a flexible method to estimate 3D space use probabilities, which might be used with different configurations of acoustic receiver arrays and could even be extended to other Eulerian-based dataset collections (e.g. camera traps, RFID reader) or coarsely sampled telemetry data (e.g. geolocators). By incorporating the topographical data, our method provides a geographically and temporally explicit view of the space use of animals moving in complex 3D environments. In future, comprehensive movement models (continuoustime, state-space) might also be used for the generation of synthetic paths, adding more mechanistic insights into the resulting outputs.

ACK N OWLED G EM ENTS
We are grateful to X. Roijals for his technical support at the

DATA ACCE SS I B I LIT Y
The FISHTRACK3D package, all the documentation and a sample dataset are hosted at https ://github.com/aspil laga/fisht rack3 d/ F I G U R E 5 Overlap indices between volumetric (3D-UD, a-b) and 2dimensional (2D-dBBMM, c-d) space use estimations of common dentex individuals for consecutive weeks before and during the spawning period. Panels (a) and (c) correspond to 50% UD contour levels, and panels (b)