A quantitative, hierarchical approach for detecting drift dives and tracking buoyancy changes in southern elephant seals

Foraging behaviour of marine predators inferred from the analysis of horizontal or vertical movements commonly lack quantitative information about foraging success. Several marine mammal species are known to perform dives where they passively drift in the water column, termed “drift” dives. The drift rate is determined by the animal’s buoyancy, which can be used to make inference regarding body condition. Long term dive records retrieved via satellite uplink are often summarized before transmission. This loss of resolution hampers identification of drift dives. Here, we develop a flexible, hierarchically structured approach to identify drift dives and estimate the drift rate from the summarized time-depth profiles that are increasingly available to the global research community. Based on high-resolution dive data from southern elephant seals, we classify dives as drift/non-drift and apply a summarization algorithm. We then (i) automatically generate dive groups based on inflection point ordering using a ‘Reverse’ Broken-Stick Algorithm, (ii) develop a set of threshold criteria to apply across groups, ensuring non-drift dives are most efficiently rejected, and (iii) finally implement a custom Kalman filter to retain the remaining dives that are within the seals estimated drifting time series. Validation with independent data sets shows our method retains approximately 3% of all dives, of which 88% are true drift dives. The drift rate estimates are unbiased, with the upper 95% quantile of the mean squared error between the daily averaged summarized profiles using our method (SDDR) and the observed daily averaged drift rate (ODDR) being only 0.0015. The trend of the drifting time-series match expectations for capital breeders, showing the lowest body condition commencing foraging trips and a progressive improvement as they remain at sea. Our method offers sufficient resolution to track small changes in body condition at a fine temporal scale. This approach overcomes a long-term challenge for large existing and ongoing data collections, with potential application across other drift diving species. Enabling robust identification of foraging success at sea offers a rare and valuable opportunity for monitoring marine ecosystem productivity in space and time by tracking the success of a top predator.

Foraging is a central element of an animal's life. Being a successfully forager is directly translated into survival, reproduction, and ultimately population growth 1 . Foraging activity (where, when and how individuals acquire resources), is therefore a core concern that underpins ecological research. Acquiring this information from terrestrial systems is difficult but tractable. However, collecting information on foraging behaviour from marine animals is especially challenging because their oceanic environment limits our ability to make direct observations of feeding activities.
Broadscale approaches to studying the foraging ecology of marine predators include stomach contents 2 , stable isotopes 3 , fatty acid signature 4,5 and genetic methods 6 . Animal telemetry approaches, with the on-going development and miniaturization of sensors, provides increasingly detailed insight into many aspects of marine www.nature.com/scientificreports www.nature.com/scientificreports/ post breeding and one during the post-moult trip (see [55][56][57][58] for a full description of the general field procedures).

Seal instrumentation was carried under ethics approval form the Australian Antarctic Animal Ethics Committee (AAS 2265 & AAS 2794) of the Australian Antarctic Division.
From all the dives (n = 18064) recorded by the tags, we only kept those dives reaching a minimum depth of 100 m, to avoid potential bias from residual air in the lungs. We also removed those dives with a duration lower than 300 sec as they are short, shallow, exploratory dives 39,44,59 that resulted in a final set comprising 97.5% of the original dives (n = 17622, Table 1). We visually inspected all dives meeting these criteria and classified them as potentially drift or non-drift dives according to their shape and the velocity records; potential drift dives requiring a passive ascent or descent phase, without directional changes ("wiggles") ( Fig. 1) or abrupt changes in the recorded velocity. Drift dives were further allocated as certain or uncertain, and as positive or negative. We then summarized each high-resolution VTDR dive using a Broken-stick algorithm (BSA). This reproduces the on-board processing of dives that occurs on the SRDLs 49,50 , resulting in a summarized form with only four subsurface inflection points retained together with the start and end points (Fig. 1).
Drift dive selection process. For the automated drift dive selection process on the summarized dive profiles, we introduce two new important steps. First, we pay particular attention to the order in which the inflection points are selected by developing a 'Reverse' Broken-stick algorithm (hereafter RBSA), similar to a recently implemented approach 49  www.nature.com/scientificreports www.nature.com/scientificreports/ with respect to dive profile characteristics to apply across groups, to automatically select those dives whose drift rates will be submitted to the final Kalman filtering stage.
Reverse Broken-Stick algorithm (RBSA). The RBSA generates the order in which the inflection points are selected before the satellite transmission. The basis for the RBSA is the same as the original BSA. From the summarized profiles the only inflection point whose original position is known is the deepest point, that is, the first selected point (as it shows the highest difference between the original high-resolution time-depth profile and the surface). A linear profile is constructed between this deepest point and the start and end points of the dive. The second point will then be determined by the largest discrepancy between the reconstructed profile and the transmitted points. The RBSA recursively reconstructs the profile until the fourth and final point is found 49 . The RBSA also generates the original residuals from the BSA the last of which, i.e. the largest remaining difference between the summarized and detailed profile, gives a relative indication of the amount of vertical activity not well captured by the summarized profile.
The inflection points are transmitted sorted by time of occurrence along the dive, not by the order in which they are selected by the BSA. Thus, an inflection point order of [2.1.3.4] indicates that the first point selected by the BSA will be the second timestamp (T2, D2) of the dive profile, the second inflection point selected is the first timestamp (T1, D1) and so on. This inflection point order (ifp) is used to organise the dives into groups.
The distribution of known (certain) drift dives from the high-resolution VTDR data was checked across these RBSA groupings, and used to make a first pre-selection of candidate drift dives (i.e. those dives that may be drift dives). Of the possible 24 groups identified by the inflection point ordering, eight comprised the majority of the known drift dives (>90% total). These groups became the only ones considered as potential drift dives (step 1, Fig. 2) and retained for the following calculations and selection procedures.
Drift rate estimate. A combination of different dive proportions and the position of the deepest point are then used to determine the drift segment of each potential drift dive (Appendix A and Table 2). Drift rate (m s −1 ) is then calculated as the difference of depth between the start and end point of the drift segment with respect to the time spent on the segment, i.e. Dr = ∆(D)/∆(T) (see Fig. 1). The sign of the drift rate allows us to further allocate dives into negative/positive subgroups.
Developing threshold criteria. For each individual dive profile, a number of numerical variables were calculated based on dive proportions. Proportional values were used to scale variables irrespective of shorter/longer or shallower/deeper absolute profiles, to minimize the influence of seals diving variability. Figure 1 provides a visual aid for these variables, which each give information about the dive shape, and a detailed description of the threshold criteria is provided in Appendix A. For each dive group we separately constructed density plots of these variables for both certain drift and non-drift dives, allowing different threshold criteria to be developed and automatically applied across the different groups (step 2 in Fig. 2).
The criteria selection and its thresholds for each group were developed sequentially as follows: First, density plots of all the proposed criteria for each seal were constructed, to show the degree of overlapping for drift and non drift dives. The criterion with the lowest degree of overlap was then selected first. This selected criterion was inspected in greater detail with an accepting-rejecting plot (see Appendix A) to find an optimal threshold; aiming for a reduction of ~50% of the non drift dives at a cost of losing as much as 5-10% of the true drift dives. Once the first threshold was identified, it was applied to the dataset, and the previous step was iterated for each of the 15 groups until no further optimal threshold could be found.
Kalman filtering drift rates. The two-step process described above supplies a final set of candidate drift rates to a custom Kalman filter. This is implemented to remove those dives with unrealistic drift rates in relation with the seal drift time series.
Step 3 in Fig. 2 gives a schematic representation of the filtering process. Kalman filters are a family of methods used to filter time series and reject/recalculate points using the trajectory of the signal, for example to filter noisy animal movement paths 60 . We applied the Kalman filter to the drift rate time series of each individual animal. Our Kalman filter assumes that (i) the vertical drift rate of a seal is proportional to the squared root of the difference between water density and the seal body density, (ii) water density is constant, and (iii) seal density changes through mass accretion (primarily blubber), and takes the general form: where: ρ k = seal density at dive k m 0 = initial mass δ k = mass increment v 0 = initial volume Vδ k = increment of volume associated with the mass increment at dive k μ k = buoyancy of the seal at dive k α = constant sign(ρ k − 1) = sign of the seal's density (it is lost on the square root calculation) ρ − 1 k = squared root of the difference between seal density and water density at dive k (2019) 9:8936 | https://doi.org/10.1038/s41598-019-44970-1 www.nature.com/scientificreports www.nature.com/scientificreports/ We then modelled the mass increments δ k as a random walk: where: δ k = mass increment at dive k δ k−1 = mass increment atdive k − 1 η k = variation on the mass increment associated with the dive k τ (δ) = variance of η k dependent on the masss (t k − t k−1 ) = timeincrement between the current dive and the previous one and we consider the error on the drift rate observations as normally distributed: www.nature.com/scientificreports www.nature.com/scientificreports/~μ r k = observed drift rate for the dive k μ k = buoyancy of the seal at dive k z r ( ) k τ = variance on the observed drift rate conditioned on whether or not dive k belongs to the trajectory The Kalman filter evaluates whether any drift rate observation associated with the time-varying density change process is inside or outside the most likely trajectory of the time series based on the expected variation associated with both the process and the observation. Whether any potential drift dive is inside or outside the trajectory of the drifting time series is defined as a binomial state variable Z k with two possible outcomes: 1 (dive k is inside the trajectory) or 0 (dive k is outside the trajectory):z The most probable drift rates of observations that are unlikely to be inside the trajectory of the drifting time series may also be estimated. However rather than using the estimated drift values, we accept as drift dives only those with a probability of being inside the trajectory close enough to 1 [P(Z k = 1) > 0.95] and retain these observed drift rates.
Validation of the method. Drift rate evaluation. To obtain the "true" drift rates we computed the rate of change in depth for all time-steps inside each drifting segment from the original high-resolution time-depth records. To check the robustness of our final drift rates the median of all these values was then subtracted from the value extracted from the summarized dive profile as a measurement of bias.
We also generated an observed daily averaged drift rate (ODDR) from the high-resolution drift dives. These were compared with a daily averaged drift rate calculated from the summarized profiles using our method (SDDR), and the difference between both used as another direct measure of bias.
To quantify the improvement in estimations due to the Kalman filter we calculated the SDDR before and after the application of the Kalman filter. We compared the performance by calculating the squared error and its mean (msr) between both SDDR's and the ODDR.
Validation with independent seal data. To assess our model performance we processed six additional data sets from Macquarie Island deployments during 2004 post-breeding trips as well as four from 2005 post-moulting trips (n = 10 seals). We visually inspected the high-resolution profiles of those dives accepted by our hierarchical procedures, and calculated the proportion that were true drift dives.

Results
From the 17622 high-resolution dive profiles visually classified, 1072 (6.1%) and 250 (1.4%) were classified as certain or uncertain drift dives, respectively (Table 1). Seal b14303, undertaking a post-moult trip, was the only one with identifiable certain positive drift dives, with a proportion of 4.8% and 1.6% for negative and positive drift dives respectively (Table 1). We found an average of 5.1% (range 4.6-6.7%) of certain negative drift dives across the three seals (Table 1). These numbers are consistent with previously reported values from this dataset 39 .  . This removed 3909 out of the 17622 dives, of which only 121 were considered certain drift dives i.e., an acceptably low (3.1%) overall false rejection rate (Appendix B). The eight retained groups together represented more than 90% of the certain negative drift dives and 86% of the certain positive drift dives (Appendix B). The criteria developed to automatically determine the drifting segment for each of the eight retained dive groups are shown in Table 2. As an example, for group [2.1.4.3] if the proportional dive time occupied by segment 1 is above 0.25, this comprises the drift segment; otherwise the drift segment comprises the next longest segment.
For each group, up to seven threshold criteria were applied sequentially to give a criteria-threshold combination that efficiently rejected certain non-drift dives. The specific criteria applied to each dive group and their threshold values are reported in Table 3. Positive drift dives occurred throughout all major groups, excluding [3.2.1.4] which represented only negative drift dives, and different criteria were applied between positive and negative drift dives within groups (Table 3). An example of a widely applied criterion is d1, the ratio between the depth of the first inflection point and the maximum depth, which for known drift dives was less than 0.6 to 0.8 across all groups. The exception was group [3.2.1.4] in which the drift segment is always segment 2 so this threshold (0.8) applies instead to the d2 criteria. After the application of the threshold criteria, 615 out of 17622 (i.e. 3.48%) candidate drift dives were retained for the Kalman filtering step.
Validation. Application of the Kalman filter rejected 155 (28.2%), 57 (34.1%) and 66 (37.3%) of the final candidate drift dives for the three test seals. The comparison of our final post-filter drift rates with the "true" rates calculated from high-resolution profiles showed there were no differences in the bias distribution among the three seals (F 238 = 0.3893, p = 0.6778, Fig. 3). The median bias after pooling the three bias distributions across the three seals was -0.0025 (S.D: 0.07, 95% confidence interval (CI): −0.03, 0.02).
Based on the comparison between the observed daily averaged drift rate (ODDR) and that obtained using our method for summarized profiles (SDDR) there were also no differences in the bias values across the three seals (F 238 = 2.61, p = 0.08). Once pooled the median value for the bias did not depart significantly from 0 (median: 0, S.D: 0.02, 95% CI: -0.001, 0.002; t 238 = 0.58, p = 0.56; Fig. 4A). There was no evidence for any trend in bias magnitude associated with an increase in the theoretical daily averaged drift rate (r = 0.02, t 238 = 0.33, p = 0.74; Fig. 4B).  www.nature.com/scientificreports www.nature.com/scientificreports/ The correlation between the ODDR and SDDR daily averaged drift rates (r = 0.99, p < 0.001; Fig. 4C) indicates our method was highly successful for the test seals.
After applying the Kalman filter, the SDDR time series efficiently followed the ODDR on all three test seals (Fig. 5). The Kalman filter implementation substantially reduced the mean squared error between the SDDR and the ODDR by an order of magnitude (95% upper CI before and after being 0.04 and 0.0015, Fig. 6).
The validation of our approach with 10 independent seals showed on average the percentage of retained dives being true drift dives was 87.5% (S.D: 9.35, Table 4).
The full filtering process as applied across the test (n = 3) and validation (n = 10) seals is visualised in Appendix C. All the filtering procedures have been implemented in R 61 and JAGS 62 and are freely available in the form of an R package (https://github.com/farcego/SlimmingDive).

Discussion
The occurrence of drift dives, where animals passively sink or rise in the water column, enables buoyancy changes to be determined in some marine species. Drift rate changes related to fat gain (or loss) provide a rare and valuable index of foraging success at sea. Here, we have presented a, reliable method to quantify drift rates from the summarized satellite relayed time-depth-record data widely used for migratory marine species. The process-based Kalman Filter is consistent with our understanding of the ecological processes governing the energy budgets of elephant seals i.e. a gain in fat during the two at-sea phases of the seals annual cycle, with effective results for both post-breeding and post-moult animals. We have not directly considered the effects of residual air as a potential source of bias on our estimates because (i) we don't consider shallow dives (i.e. less than 100 m depth) as potential drift dives, and (ii) elephant seals exhale before diving. Previous research 23 evaluated the potential effect of residual air present in elephant seal lungs, finding little effect on deep dives. That may not be the case for other shallower, breath-holding marine mammals, where these assumptions may be too strong. Our method overcomes a long-term challenge to robustly identify at-sea foraging success, and provides great opportunity for linkages between ecology, physiology, behaviour and environmental drivers to be further explored.
The new method provides a time series of drift rates, and the daily averaged values we obtain from the summarized dive profiles show good concordance with those obtained from visually inspected high resolution dive profiles. Compared with the existing approach 41 , which reported 71.4% of retained dives as being true drift dives, our approach retained 87.5%. That gives over a 15% increase in the true drift rate retention, reducing the impact of false positives on the estimated drift rate time series, and contributing to reduce the error/variance of the estimation.
Inclusion of false positive drift dives can result in a higher variance among drift rate estimates and require some further processing; often achieved with smoothing/interpolating techniques such as splines 23,41,59 . Such smoothing/interpolating techniques are based on purely statistical approaches, without any biological process underpinning them. Using a custom Kalman filter incorporates a biologically relevant mechanistic model. Although this filter does not remove every non-drift dive, it greatly reduces their occurrence to approximately 10% of the final set of retained dives. The filter also reduces the variability of the daily drift rate estimates, by over   An important improvement from previous approaches is that our method can detect when the seal is positively buoyant. Positive buoyancy has implications for quantifying the individual foraging behaviour and success of individuals, as well as the quality of the foraging grounds. In our study, we have processed five post moulting trips, of which three exhibited substantial periods of positive buoyancy up to 150 days, Compiling a realistic record of daily body condition changes would have not been possible with previous approaches which is the ultimate goal of our approach. We make our method available to the research community in the form of an R package under a General Public License.
The results are also consistent with expectations regarding the energy budgets of seals 23,44,59 . All the seals exhibit their lowest body condition at the start of the foraging trip, after fasting for 1-2 months (Appendix C). They show a progressive increase of body condition as they remain at sea, indicating that they are foraging sufficiently well for their physiological needs and to gradually replenish their lipid reserves. In the longer post-moult time-series periods of positive buoyancy occur, often followed by a return to negative buoyancy.
Once buoyancy changes at temporal scales of days to months for individual animals can be estimated these data can be used to relate patterns of individual foraging success to factors such as such as who lives or dies, or who pups successfully, and how this links to where (spatially) and how (functionally) individuals may forage. Compiling patterns of foraging success across individuals will facilitate population level studies such as why some populations are stable and others declining [63][64][65] . Southern elephant seals have been tagged from all Southern Ocean breeding populations 35 , a global effort spanning more than two decades. Many hundreds of individual animals have been tagged, including both sexes as well as adults and juveniles 66 . Our automated approach is tractable for analysing existing and ongoing large dataset collections for larger overarching studies for example that link   Table 4. Validation of the drift dive methodology with 10 independent Macquarie Island seals. Seal id = reference code for each individual tag/seal. Trip: pb = post-breeding trip, pm = post-moulting trip. N = Total number of dives recorded by each tag. Rd = number of retained dives after the application of our method. %d = proportion of dives retained from the total number of dives recorded. Rdd = number of retained drift dives. %Dd = proportion of the retained dives that were true drift dives, as determined by visual inspection of all retained dives using the original high resolution time-depth profiles.