nent example of the latter outcome (D'Agrosa et al. 2000, Jaramillo-Legorreta et al. 2007, 2019).

Short-beaked common dolphins *Delphinus delphis* Linné, 1758 in the Bay of Biscay (northeast Atlantic) have been subjected to potentially dangerous levels of bycatch since the 1990s (see Murphy et al. 2019 for a review). In France, the magnitude of bycatch is evidenced by high numbers of strandings of bycaught dolphin carcasses on the Atlantic seaboard. Mannocci et al. (2012) estimated that with 1000 bycaught individuals per year, the northeast Atlantic population could be reduced by 20% in 30 yr and could become extinct in a hundred. Still, the phenomenon has become more intense in the past decade. Since 2016, the number of strandings has increased every year (http://seamap.env.duke.edu). Bycatch estimates from strandings that occurred solely in the Bay of Biscay have been above what the population of common dolphins of the whole northeast Atlantic can theoretically sustain (Mannocci et al. 2012, Peltier et al. 2016, ICES 2020). From January to April 2019 alone, strandings suggested bycatch levels as high as 9500 individuals (95%CI: [6890; 14 200]) (Dars et al. 2019).

Stranding data can be complementary to data from European observer programs and can provide relevant information on cetacean bycatch (Peltier et al. 2016, IJsseldijk et al. 2020). In France, they have been instrumental to estimate the magnitude of common dolphin bycatch in the Bay of Biscay (Peltier et al. 2016) and to reveal a strong seasonality, with increased mortality in winter months (January to March) (Peltier et al. 2016, 2019, 2020, Dars et al. 2019). From 2016 onward, strandings of animals showing bycatch evidence have risen in summer (July and August), although to a lesser extent than in winter (Dars et al. 2019). Stranding data further made it possible to identify spatiotemporally varying mortality areas (as areas where carcasses of bycaught dolphins were released from fishing boats) through reverse drift modelling of carcasses trajectories (Peltier & Ridoux 2015, Peltier et al. 2016, 2021). In winter months, most bycatch events likely occur in the southern part of the continental shelf and slope (Peltier et al. 2016, 2020), resulting in strandings all along the shore from the south of the Finistère region to the border to Spain, whereas in summer most carcasses strand in southern Finistère (Fig. 1).

Finally, the bycatch pattern displays high interannual variability, with different stranding levels and associated bycatch level estimates from year to year (Peltier et al. 2016).

Stranding data further revealed possible association between bycatch events and common dolphin preys. The reverse drift modelling of carcasses allowed the identification of some potentially high-risk fisheries and their caught fish species (Peltier et al. 2020, 2021). The analysis of stomach content of bycaught stranded animals suggested that most dolphins were feeding when death occurred. Moreover, some fish species targeted by common dolphins (e.g. mackerel or sardine) are also included in the diet of predatory fish species of fisheries interest (e.g. hake or seabass) (Spitz et al. 2013). Hence, changes in the local distribution and abundance of prey species might be a substantial driver of the co-occurrence of common dolphins, commercial fish species and fisheries in localized areas.

As marine mobile predators, cetaceans have a dynamic distribution integrating ecological processes across all levels of the trophic web (Croll et al. 1998, Barlow et al. 2020). The proximal relationship between prey and predator could drive the spatiotemporal variability of the bycatch pattern: fish distribution is notoriously variable in both space and time (Hyrenbach et al. 2000) and is governed by many factors ranging from dynamic oceanographic conditions to each species' annual cycle. In this study, we aimed at assessing the influence of oceanographic

<figure><img src="image_159.png" /><figcaption>Fig. 1. Study area showing bathymetry and geographic references mentioned in the text, and location of common dolphin stranded carcasses used for the study (red dots)</figcaption></figure>
processes on bycatch of common dolphins in fishing gear. We thus focused on the distal relationship between marine predators and oceanographic processes. We assumed oceanographic processes spatiotemporally structure the availability and aggregation of prey, creating areas prone to attract common dolphins and predatory fish targeted by fisheries and consequently increase the bycatch risk. For instance, meso and submesoscale processes such as fronts (where 2 joined water masses differ by their density), upwelling or eddies enhance the enrichment, concentration and retention of nutrients (Bakun 1996). These processes facilitate the development of trophic networks, from phytoplankton to zooplankton, to fish and finally, to apex predators (Bakun 1996, 2006).

Knowledge of oceanographic conditions is readily available from either remote sensing or modelling. This information is therefore commonly used to make inferences about prey availability and species distribution (Forney 2000, Becker et al. 2010, Best et al. 2012, Stephenson et al. 2020). Oceanographic data have been previously used in the context of bycatch studies, either to identify high risk areas through the association of species distribution modelling and distribution of fishing effort (Žydelis et al. 2011, Murray & Orphanides 2013, Di Tullio et al. 2015, Díaz López et al. 2019) or to identify possible drivers of bycatch events (Gardner et al. 2008, Cosandey-Godin et al. 2015, Hahlbeck et al. 2017, Scales et al. 2018). Cosandey-Godin et al. (2015) suggested the forecasting of high bycatch risk areas to identify mitigation measures. We investigated relations between mortality of common dolphins in fishing gear and oceanographic processes in the Bay of Biscay. Predicting high bycatch areas was beyond the scope of this paper, although it constitutes a promising venue for further investigations.

This work capitalized on 2 sources of data: the bycatch-mortality index inferred from strandings and fine-scaled essential oceanographic variables (EOVs) (García-Barón et al. 2020, Tew-Kai et al. 2020). We selected physical variables significantly contributing to the characterization of the ocean realm (Tew-Kai et al. 2020). We used hierarchical Bayesian models to account for spatiotemporal processes: bycatch mortality areas, oceanographic processes, and common dolphin and fisherie distributions are dynamic in both space and time. The area of study was limited to waters of the Bay of Biscay, with a small part of the Celtic Sea (Fig. 1). After model fitting, we used out-of-sample, 1-step-ahead prediction to predict bycatch mortality in 2019 from oceanographic processes.

## 2. MATERIALS AND METHODS

### 2.1. Dolphin mortality index from stranding data

#### 2.1.1. Mortality index

For most bycatch events, no data or location at sea is available. However, some bycatch events are detectable from the stranding of common dolphins presenting bycatch evidence (Kuiken 1994). To determine the origin of carcasses, a reverse drift trajectory can be computed with the drift prediction model MOTHY (Modèle Océanique de Transport d'Hydrocarbure, developed by Météo France; Daniel et al. 2002), integrating drift conditions and characteristics of dolphins' carcasses (see Peltier & Ridoux 2015, Peltier et al. 2012). Reverse drift trajectories of bycaught animals were used to compute a mortality index (MI) to identify probable mortality areas and to quantify the intensity of mortality events at sea in space and time (Fig. 2).

From 2012 to 2019 stranding data, only carcasses presenting bycatch evidence or which stranded during an intense mortality event were selected. Intense mortality events are when more than 30 dolphins strand on 200 km of coastline within 10 d (Peltier et al. 2016). We further selected carcasses with estimated mortality windows up to 15 d (2384 out of 4544, i.e. 52%). For each carcass, a reverse drift trajectory was calculated with MOTHY, with one point every 10 h. Drift tracks were limited from 0 to 5 or from

<figure><img src="image_160.png" /><figcaption>Fig. 2. Spatialized mortality index (MI) for March 2018. MI was computed from reverse drift modelling of stranded carcasses of common dolphins showing bycatch evidence. It allows us to locate potential mortality areas and to quantify the intensity of mortality events</figcaption></figure>
5 to 15 d depending on the mortality window estimated from the decomposition state of the stranded carcasses. The computation of the MI on a daily basis, using a regular grid of 0.1° × 0.1° resolution (i.e. approximately 26 nmiles² or 89 km²) for the Bay of Biscay, assumed an equal probability of bycatch occurring at each drift point. For a day d, a grid cell c and a carcass i, MI was defined as follows:

$$MI_{d,c,i} = \frac{\text{No. of drift points of } i \text{ in grid cell } c}{\text{Total no. of drift points of } i} \quad (1)$$

The more a carcass stagnates in a grid cell, the higher MI reflects a higher probability that death occurred in that grid cell. For a given cell on a given day, there could be reverse drift points associated with several carcasses. The total MI for d and c was then:

$$MI_{d,c} = \sum_{i} MI_{d,c,i} \quad (2)$$

The greater the number of carcasses that have drifted through a grid cell, the higher its MI. Finally, daily MIs were summed for each grid cell over each month from 2012 to 2019 resulting in the following spatial monthly MI, for month m and grid cell c (Fig. 2):

$$MI_{c,m} = \sum_{d \in m} MI_{d,c} = \sum_{d \in m} \sum_{i} MI_{d,c,i} \quad (3)$$

The total sum of MI over each month is equal to the number of stranded dolphins included in the dataset. In total, 1789 carcasses stranded all along the shore of the Bay of Biscay (Fig. 1) were used for the present study, with substantial seasonal and inter-annual variations (Fig. 3). Monthly numbers are provided in Supplement 1; www.int-res.com/articles/suppl/m679 p195_suppl/.

### 2.1.2. Stranding probability

A wide range of factors lead a carcass to strand, including buoyancy (the carcass must be floating) and drift conditions (wind being one of the main drivers). Consequently, not all carcasses of bycaught common dolphins do strand (Peltier et al. 2016). Drift conditions can sometimes keep carcasses at sea rather than push them to the coast. Buoyancy of carcasses is probably influenced by several parameters (i.e. depth at which death occurred, buoyancy at time of death; Moore et al. 2020) so that they may either sink or float when they are released from fishing boats. The calculated MI thus contained ‘false’ zeros, associated to grid cells where conditions would not have led to strandings if carcasses happened to drift through

them. Conversely, ‘true’ zeros were associated to grid cells where conditions could have led drifting carcasses to strand. To account for this, we systematically included stranding probability as a covariate in all models. Stranding probabilities were estimated from drifts of simulated dolphin carcasses with MOTHY (Peltier et al. 2013, 2014, 2016). It was calculated over 10 d periods and averaged over each month. A conditional autoregressive (CAR) model was then applied to the result to leverage spatial dependence and smooth estimated probabilities. The spatial resolution of estimated probabilities was larger than that of MI (Fig. 4): each grid cell of the MI grid was matched to its closer neighbour in the stranding probability grid (using R package sf of Pebesma 2018). True zeros should be associated with high stranding probabilities whereas false zeros should be associated with low stranding probabilities.

## 2.2. Oceanographic data

Monthly EOVs and derived variables were computed from the 3-dimensional HYCOM numerical model of the French Service Hydrographique et Océanographique de la Marine (Shom). Based on the original HYCOM model (Bleck 2002), it was set with a number of numerical developments to optimize the model for coastal zones (Pichon & Correard 2006, Morel et al. 2008, Lahaye et al. 2011) and has already been successfully used over the Bay of Biscay area (e.g. Pichon et al. 2013). Hindcasts were calculated hourly from 1 January 2012 to 31 December 2019 to compute monthly EOVs and derived variables. EOVs are outputs of the model (temperature, salinity, cur-

<figure><img src="image_161.png" /><figcaption>Fig. 3. Evolution of the total mortality index (MI) per month and per year. The monthly sum of spatial MI over the study area is equal to the number of stranded dolphin carcasses for which a reverse drift was calculated, thus giving a mortality signal</figcaption></figure>
Author copy

<figure><img src="image_162.png" /><figcaption>Fig. 4. Stranding probability in March 2018, estimated from simulated drifts of common dolphins, then smoothed with conditional autoregressive (CAR) model to account for spatial dependency between adjacent cells</figcaption></figure>

rent), and derived variables have been calculated from these EOVs. One EOV and 2 derived variables, hereafter referred as oceanographic variables (OVs), were integrated in this study: (1) the sea surface temperature (sst), (2) the mean sea surface temperature gradient (mean_sstGrad), quantifying the intensity and wideness of thermal fronts, and (3) eddy kinetic energy (eke), quantifying the turbulent component of the residual (non-tidal) surface current. sst was chosen for its strong seasonal component (Lambert et al. 2017) but also as it is commonly integrated in habitat modelling of multiple marine taxa (Doniol-Valcroze et al. 2007, Mellin et al. 2012, Hughes et al. 2014, Castro et al. 2020). In addition, thermal fronts and eddies have been identified as key structures for the conservation of marine predators (Scales et al. 2014, 2018). Such interfaces might trigger bottom-up processes. As nutrient mixing and retention in the photic zone can boost primary productivity, interfaces like fronts and eddies might attract apex predators through a cascading effect from the lower trophic levels to the highest ones (Bakun 1996, 2006, Scales et al. 2014). This association was observed for several marine predators (Nel et al. 2001, Doniol-Valcroze et al. 2007, Bailey & Thompson, 2010, Miller et al. 2015, Snyder et al. 2017) along with further associations with bycatch events (Hahlbeck et al. 2017, Scales et al. 2018). Moreover, the Bay of Biscay encompasses a

shelf area with complex dynamic mesoscale processes following a seasonal cycle (Yelekçi et al. 2017), and we meant to integrate some of these features in our analysis through the dynamic variables eke and mean_sstGrad.

As the spatial resolution of the oceanographic data was the nautical mile, it was upscaled to match that of the MI. All data can be visualized through maps and histograms at http://pelabox.univ-lr.fr:3838/pelagis/DdSeaByc/ (Tab: '1. Data').

## 2.3. Spatiotemporal modelling framework

Our goal was to elucidate if oceanographic conditions were associated with bycatch of common dolphins in fishing gear with statistical modelling. We considered a general set up MI ~ OVs where MI is the response variable whose variations are to be correlated with OVs, with linear effects for covariates.

### 2.3.1. Spatial and temporal considerations

Common dolphin bycatch—and stranding of animals showing bycatch evidence—is a seasonal phenomenon in the Bay of Biscay (Dars et al. 2019, ICES 2020). The mortality signal, defined here as the number of stranded common dolphin showing bycatch evidence, is higher in January to March (winter) and in July and August (summer) but is lower in between (Fig. 3). Oceanographic conditions also follow a seasonal cycle with different dynamics in winter and summer and transition regimes in between (Yelekçi 2017). Our exploratory analysis (not shown here) suggested that the relationship between the response variable and oceanographic covariates varies across months. Hence, we chose to fit one model per month to get a climatology over the study period. Each monthly model gathered data for the corresponding month from 2012 to 2019, resulting in 12 models (from January to December). A similar approach was used in Foravanti et al. (2021). This climatological approach had the advantage of reducing computational cost.

The relation between OVs and MI in one given month can vary from one year to the next, for in-
stance because of a timelag in cyclic phenomena (Huret et al. 2018) or because of singular climatic events. The dynamic of mesoscale processes is complex in the Bay of Biscay (Yelekçi 2017). Singular phenomena could be especially attractive for preys of dolphins and predatory fish targeted by fisheries depending on their extent, their duration and their time of occurrence. The effect of OVs on bycatch risk (and thus on MI) must be varying in time, as suggested by the intensification in strandings from 2016 (Dars et al. 2019). This led us to consider annually varying coefficients for OVs (i.e. random slopes). We further added an annually varying intercept (random intercept) to account for different annual levels of bycatch. Using random slopes amounts to considering an interaction with year: time-series of the effect of OVs on MI for each month and each year were an output from models over the study period.

Generalized linear models (GLM) and generalized additive models (GAM) usually assume all relevant covariates have been included in the model (Zuur et al. 2007a,b). This assumption is unrealistic in most cases, and omitted variables may result in a residual dependence that, if not accounted for, may result in misleading inferences (Legendre et al. 2002, Valcu & Kempenaers 2010). Omitted variables bias may manifest in residual autocorrelation, either temporal, spatial, or both. Because there was one model per month, the time lag between each 'measurement' was 11 mo, and we assumed no temporal autocorrelation. Yet, spatial autocorrelation, if not taken into account, could suggest spurious correlations with OVs. Consequently, all models included explicitly a spatial component in addition to the (implicitly spatial) OVs.

## 2.3.2. Hierarchical Bayesian framework and Gaussian random field

We used a hierarchical Bayesian framework as they are adapted to the modelling of complex spatial or spatiotemporal phenomena (Cosandey-Godin et al. 2015, Martínez-Minaya et al. 2018). Including a spatial random effect, e.g. with a Gaussian random field (GRF), is straightforward (Blangiardo et al. 2013, Cosandey-Godin et al. 2015). The GRF, indexed on each grid cell, is a stochastic process that will gather all the spatially implicit information in the data that are not accounted for by other model components, in our case oceanographic covariates (Cosandey-Godin et al. 2015). To model this random field, we chose a

CAR process (Besag et al. 1991, Besag & Kooperberg 1995) based on the notion of nearest neighbours (the so-called first law of geography: 'Everything is related to everything else, but near things are more related than distant things'; Tobler 1970, p. 236). For the vector of grid cells ($η_1, η_2, ..., η_n$) here n = 3487) it can be written as follows:

$$η_i | η_j, i \neq j, \tau \sim N\left(\frac{1}{n_i}, \frac{1}{n_i \times \tau}\right) \quad (4)$$

where $n_i$ is the number of neighbours for grid cell i and $τ$ is a precision parameter. i ~ j means that the 2 grid cells i and j are neighbours.

## 2.3.3. Model set up

MI was first transformed to $log(MI + 1)$ and then modelled by a Gaussian distribution with an identity link function. The spatial CAR process was considered constant over years, so that there would be one spatial field per month. The stranding probability (SP) was always included as a covariate with a linear effect. For oceanographic variables effects and for the intercept, we included the possibility for the coefficients to vary from one year to the next by specifying a Gaussian random (independent and identically distributed) effect on years (random slopes). For an overview of our monthly model, we had for the response variable MI in year t and grid cell c, with oceanographic covariates X:

$$ \left\{ \begin{array}{l} \text{For all } m \text{ in } 1:12 \text{ (dropping the } m \text{ subscript below for clarity)} \\ \log(\text{MI}_{c,t} + 1) \sim N(\pi_{c,t}, \sigma) \\ \text{Id}(\pi_{c,t}) = \gamma \times \text{SP}_{c,t} + \beta_{0,t} + \sum_j \beta_{j,t} \times X_{j,c,t} + \eta_c + \varepsilon_{c,t} \end{array} \right. \quad (5) $$

where N defines the normal distribution of location parameter $π_{c,t}$ and scale parameter $σ$, Id is the identity link function, $γ$ the coefficient associated with SP, $β_{0,t}$ the intercept for year t, $β_{j,t}$ the annual coefficient of covariate $X_{j,c,t}$ at cell c in year t and $η_c$ the spatial effect associated to grid cell c, which is assumed to result from a CAR model. Finally, $π_{c,t}$ are unstructured errors, i.e. residuals. Parameters estimated in the linear predictor, i.e. $γ, β_{0,t}, β_{j,t}, η_c$ and $π_{c,t}$ constitute the latent field. They are assumed to follow Gaussian distributions (Rue et al. 2009). These parameters depend on a vector of hyperparameters $θ$ that are not necessarily Gaussian (Rue et al. 2009). In this manner, the monthly model defined in (Eq. 5) is a latent Gaussian model (LGM) (Rue et al. 2009). It presents a spatial component ($η_c$) along with annual components ($β_{0,t}$ and $β_{j,t}$).
### 2.3.4. Inference with integrated nested Laplace approximations (INLA)

Inference for such complex spatial models is computationally demanding and requires approximation to be fully scalable with a large amount of data. For a large dataset with 3487 grid cells observed 8 times (i.e. years) per analysis, i.e. 27 896 observations per month, classical methods like Markov Chain Monte Carlo may be prohibitively slow (Blangiardo et al. 2013). We used INLA, which offers a very efficient alternative giving fast and reasonable estimations (Rue et al. 2009). Especially adapted to LGMs, this algorithm uses a combination of analytical approximations and numerical integration for the inference of posterior marginals and likelihood, which greatly reduce computational demand while being extremely accurate (Rue et al. 2009, Moraga 2019). Inference with INLA was carried out with the R-INLA package in R statistical software v4.0.2 (R Core Team 2020).

Before model fitting, all oceanographic variables were standardized (mean-centered and scaled to unity), and eke and mean_sst_grad were log-transformed to account for skewness and strict positivity. Weakly informative priors were chosen so that estimations would mostly be driven by the data (Gabry et al. 2019). Fitted and predicted values that were negative after back transformation were set to zero.

### 2.3.5. Model selection

All possible combinations of the different effects and covariates were tested, with the SP systematically included as a confounding factor (list of tested models is provided in Supplement 2). The first combinations were the simplest models with only linear effects for oceanographic covariates (i.e. in Eq. 5, $\beta_j$ replaces $\beta_{j,t}$) and no spatial field (8 possibilities with the null model). Second, coefficients were allowed to vary between years for oceanographic covariates (7 possibilities) with random slopes. Finally, the CAR process was added, with and without oceanographic covariates, considering strict linear effects or annual random slopes (15 possibilities). The 30th model was the more complex, defined in Eq. (5), where $j \in [1,3]$ and $X_j$ are thus sst, eke or mean_sst_grad. Model selection was conducted using the Wanatabe Akaike information criterion (WAIC) (Watanabe 2010). WAIC

gives an approximation of the out-of-sample predictive accuracy using the in-sample data (Gelman et al. 2014).

### 2.3.6. Model checking

The coefficient of determination $R^2$ for each monthly model was computed. $R^2$ is the proportion of variance in the response variable accounted for by explanatory variables (Johnson 2014). We appropriated the total variance to different model components, i.e. spatial or annual, and assessed their respective contributions. For random slopes models, some variance components depend on observations (Johnson 2014). For instance, the variance of $\beta_{sst,t}$ coefficient of variable sst on year $t$, depends on values of sst during year $t$. We used the matrix formula (11) of Johnson (2014); see Supplement 3 for further details.

Finally, we assessed the predictive power of models, first by conducting out-of-sample cross-validation (CV), then by testing different prediction scenarios. While the dataset was complete from 2012 to 2018, for 2019 only values of oceanographic variables were given to get predictions of MI for this year. As the intercept and covariates coefficients were estimated annually, there was an index associated to each year. For the CV procedure, we considered 2019 as a new year (predicted from the posterior distribution) and attributed to this year a new and unique index. In addition, we tested different repetition scenarios (RS): instead of attributing a new index to year 2019 and consider it as a new sample, we considered 2019 as a repetition of a previous year by attributing to 2019 its index. This way, all model components fitted during the 'repeated' year were used (covariates' coefficients, effect of the stranding probability, spatial CAR process) for the prediction of MI in 2019, still considering 2019's oceanographic conditions. We tested all 7 possibilities, i.e. we considered 2019 as a repetition of 2012, then of 2013, 2014, 2015, 2016, 2017 and 2018. Then, we confronted the predicted total mortality for 2019 to 2019's total mortality from the dataset (never seen by the models). We also computed a root mean squared error (RMSE) per RS and for the CV, based on the monthly differences between predicted and observed (from strandings) total monthly mortality. With these repetition scenarios, we aimed at addressing whether oceanographic conditions of past years affected bycatch mortality risk in the same way and could predict the observed mortality in 2019.
# 3. RESULTS

## 3.1. Model selection

For most months, the most complex model detailed in formula (Eq. 5) had the smallest WAIC values (table with WAIC of all tested models provided in Supplement 4). It included a CAR process, year-varying intercepts and year-varying slopes for the 3 covariates *sst*, *eke* and *mean_sstGrad*. In May, June and November the model with the smallest WAIC values excluded *eke* and for September *sst* was excluded. However, the difference in WAIC between these models and the most complex one selected for other months was small (relative difference of 0.01, 0.03, 0.1, 0.06% in May, June, September and November, respectively). Hence, we chose to keep the most complex model with all 3 covariates for all months to facilitate comparison of results over the 12 months.

## 3.2. Oceanographic covariates effects on the MI

Slope estimates of the effect of *sst*, *eke* and *mean_sstGrad* on spatiotemporal MIs were precise with tight credible intervals (see http://pelabox.univ-lr.fr:3838/pelagis/DdSeaByc/; Tab: '2. Results 2.1.'). There were no estimates for year 2019 as this year was not used in model fitting but in out-of-sample,

1-step-ahead cross-validation and predictions. Overall, effects were more intense in the winter season, from January to March (Fig. 5). They were close to zero for the rest of the year, except in July and August for 2016, 2017 and 2018 for *eke* and *mean_sstGrad*. The winter variation pattern of the effect of *eke* and *mean_sstGrad* was different every year. Their effect on MI was sometimes positive and sometimes negative, depending on years. In July and August, however, the effect of these 2 covariates was consistent from 2016 to 2018. It was negative for *eke* and positive for *mean_sstGrad*. The winter variation pattern of the effect of *sst* was consistent between years, with interannual differences in its intensity. It switched sign from one month to the next, for most years. In summer, the effect of *sst* was close to zero, except in 2017, when it was slightly negative. For all 3 covariates, 2017 emerged as a singular year.

The intercept was estimated as a mean level of MI with OVs being set at their mean values (given that they were standardized, that is mean-centered). It was greater in winter, with mostly positive values in January and February but negative values in March, except in 2017, when it was positive for all 3 months (Fig. 5). Its pattern of variation was more or less consistent over the study period, with differences in its intensity. In summer, the intercept was slightly positive only in 2016 and 2017.

<figure><img src="image_163.png" /><figcaption>Fig. 5. Means estimates from the 12 separated monthly models of intercepts and coefficients of oceanographic covariates, for each month and each year. Spatiotemporal monthly models were inferred with integrated nested Laplace approximations (INLA). *eke*: eddy kinetic energy; *mean_sstGrad*: mean sea surface temperature; *sst*: sea surface temperature</figcaption></figure>
### 3.3. Variance partitioning

The computation of $R^2$ allowed us to differentiate the contribution of the different model components to integrate MI variations (Fig. 6). The annual component included annual coefficients of covariates (random slopes) along with annual intercepts. The spatial component was the CAR process whose contribution in variance showed little between-month variation. The residual component included what was not accounted for by either of the 2 precedent components. In total, models better accounted for MI variations from January to April, the annual component being the most contributive component on this period. From January to March, 21 to 25% of the variance were accounted for by annual and spatial components. A rise in the annual component contribution in July and in the spatial component contribution in August resulted in a rise in the overall variance accounted for these 2 months. Afterward, most of the variations of MI were integrated in residuals. Note that, using formula (11) of Johnson (2014), covariances between components were included: some negative covariances terms explain that contributions of the 3 components naively sum to more than 100% for some months.

### 3.4. CV and predictions of 2019 mortality with repetition scenarios

For CV, 1-step-ahead prediction was performed. Each monthly model was calibrated on 2012–2018 data, so MI data from 2019 was not used in model fitting. 2019 oceanographic data were thus used for prediction, taking between-year variations into account. The total monthly MI of the dataset (observed) for 2019 fell within the confidence interval of total monthly predicted values, for both the winter and the summer mortality periods (Fig. 7). Total mean predicted values were below observed ones. In addition, 2019 was assumed as a repetition of each year between 2012 and 2018 (RS), to highlight if past oceanographic influences could predict observed mortality in 2019. Among the 7 RS, observed values of total monthly mortality of 2019 during the winter mortality period was within (or close to) confidence intervals of predicted values for 2 RS, when 2019 was considered as a repetition of 2013 (RS

2019 = 2013) and as a repetition of 2017 (RS 2019 = 2017) (Fig. 7). For RS of other years, predicted values and confidence intervals were lower than observed monthly mortality, from January to March. The summer mortality season was not reproduced by the models for RS of years before 2016. Overall, 2017 RS showed the lowest RMSE and gave the best prediction of 2019 monthly mortality levels (Table 1). RMSE for CV was below RMSE of RS of years 2012 and 2015 and close to that of 2018.

Additional results such as maps of fitted or predicted values and maps of the estimated spatial effects are available at http://pelabox.univ-lr.fr:3838/pelagis/DdSeaByc/ (Tab: '2. Results') for all 12 months. A temporal series of total monthly mortality of fitted values and observed values are provided in Fig. A1 in the Appendix.

## 4. DISCUSSION

We assumed bycatch risk could be correlated to the proximal relationship between the dynamic distribution of common dolphins, commercial predatory fish species targeted by fisheries and that of their common prey species, the latter being partially driven by oceanographic processes. Spatiotemporal information on at-sea bycatch mortality, fisheries, dolphins and fish prey species distribution is not easily accessible to study this proximal relationship. Thus, we used proxies, a mortality index inferred from strandings and oceanographic descriptors, to study the dis-

<figure><img src="image_164.png" /><figcaption>Fig. 6. For each of the 12 spatiotemporal monthly models inferred with integrated nested Laplace approximations (INLA), proportion (transformed with a square root) of variance accounted for by each model component. Residual: residual component; Spatial: spatial component (defined by a Gaussian random field); Annual: annual intercepts and coefficients of covariates</figcaption></figure>
tal relationship between common dolphin bycatch and oceanographic processes. We mobilized 8 yr of stranding data (from 2012 to 2019) and explored spa- tial and temporal links between basal and apex nodes of the trophic network, from oceanographic processes to short-beaked common dolphins (and fisheries). We fitted a series of 12 spatiotemporal monthly models that included 3 oceanographic vari- ables as covariates potentially contributing to the spatiotemporal structuring of bycatch risk of common dolphins in the Bay of Biscay. Results provided first insights into how the intensity of common dolphin bycatch in the Bay of Biscay might be modulated by key seasonal and dynamic oceanographic features.

## 4.1. Two seasons of mortality, two different underlying processes

Strandings highlight 2 seasons of high bycatch mor- tality in the Bay of Biscay: a main one from January to March and another (only observed from 2016 onward) in July and August (Fig. 3). Results reflected different patterns for these 2 mortality seasons, suggesting dif- ferences between the underlying phenomena leading to bycatch mortality. Outside these seasons, the ex- planatory power of oceanographic processes was very limited, and the effect of OVs on the mortality index was close to zero. Thus, it seems that OVs seasonally affect bycatch risk and should not be considered as season-invariant drivers of bycatch mortality.

### 4.1.1. Winter high mortality season, from January to March

The annual component, including OVs, accounted for most variations in bycatch risk during the winter mortality season (Fig. 6), suggesting a substantial in- fluence of the oceanographic processes. Moreover, January, February and March were the months when OVs had the strongest effects on bycatch risk (Fig. 5). Depending on OVs, their effects (i.e. monthly estimated linear coefficients) showed different inter-annual pat- terns. The month-to-month effect of sea surface tem- perature was the only one to be consistent from one year to the next. Sea surface temperature partially drives the distribution and seasonal migrations of mul- tiple marine taxa, from small pelagic fish (Lanz et al. 2009, Hughes et al. 2014,) to bigger fish (Hahlbeck et al. 2017) and high-level predators (Weltz et al. 2013). In the Bay of Biscay, the influence of sea surface tem- perature affects the distribution of identified prey

species of common dolphins and commercial preda- tory fish species (Spitz et al. 2013), such as anchovies and sardines (Planque et al. 2007, Petitgas et al. 2014, Politikos et al. 2015). Therefore, a correlation between mortality and sea surface temperature can be ex- pected if the underlying mechanisms are linked to preys of both common dolphins and commercial fish species targeted by fisheries. The effects of the mean temperature gradient and eddy kinetic energy were of the same magnitude as that of sea surface tempera- ture but presented high inter-annual variations. Mean temperature gradient quantifies the intensity and wideness of thermal fronts, whereas eddy kinetic en- ergy is associated with turbulent surface features such as eddies. Because they were averaged over each month, they should refer to relatively persistent structures. In winter, the Bay of Biscay environment is characterized by a seasonal cross-shore (west to east) surface temperature gradient with lowest temperature close to shore and intense frontal activity parallel to the coast (north to south) (Yelekçi et al. 2017). These frontal structures are freshwater fronts, correlated with the mixing of oceanic waters and cold freshwater inputs from river plumes (Yelekçi et al. 2017). Their location and timeline concur with observed patterns regarding bycatch mortality: during winter months, strandings were from all the Atlantic coast (http:// seamap.env.duke.edu), and high mortality areas were previously identified on the continental shelf (Peltier et al. 2021). Moreover, most of highest MI values were located offshore of the coast from the south of the Fin- istère region to the south French border (not shown, see http://pelabox.univ-lr.fr:3838/pelagis/DdSeaByc/; Tab: '1. Data'), especially in 2017 when strandings were the more numerous over the 2012 to 2018 period. Together with the effects of the associated OVs on the MI, results suggested these seasonal fronts may be targeted by both fisheries and common dolphins as areas where fish aggregate, thereby putting the latter at risk of bycatch by the former.

### 4.1.2. Summer high mortality season in July and August

July and August were the only 2 other months when OVs accounted for a non-negligible fraction of varia- tions in bycatch mortality (Fig. 6). The amount of total explained variance accounted for by OVs (included in the annual component) was more limited compared to the winter mortality season and was below the contri- bution of the CAR process. This could suggest that, during this summer season, other drivers (currently
omitted in the model) might drive mortality risk. How-
ever, this was expected as strandings are at least one
order of magnitude lower in summer and only be-
came substantial from 2016 onward. Nevertheless,
different patterns emerged regarding the effects of
oceanographic covariates. Conversely to the winter
season, the sea surface temperature coefficient was
close to zero for all years, and the effects of mean tem-
perature gradient and eddy kinetic energy on the MI
were, respectively, quite consistent for all 3 years when
a summer mortality peak was observed (2016, 2017,
2018). For these 3 years and in July and August, the
effect of the mean temperature gradient was slightly
positive, meaning that thermal fronts were associated
with higher MI and that of eddy kinetic energy was
slightly negative, meaning that bycatch mortality
probably occurred in low turbulence waters. In July
and August, the mesoscale dynamic activity of the
Bay of Biscay is rather different than in winter. In win-
ter, the frontal activity in the Bay of Biscay is domi-
nated by freshwater fronts resulting from input of cold
freshwater from rivers whereas in summer, there are
mainly fronts due to tidal flow (Yelekçi 2017). Winter
freshwater fronts display more spatiotemporal vari-
ability (on a given year and interannually) as they are
directly correlated to river flows (Yelekçi et al. 2017).
Summer tidal fronts are conversely quite consistent
from one year to the next because they are correlated
to a repetitive process (i.e. tides) (Yelekçi et al. 2017).
During summer, the main frontal activity is a seasonal
tidal front, called the Ushant Front, and it is located in
front of the French Finistère region (Yelekçi et al.
2017). Its activity peaks in July and August (Yelekçi et
al. 2017). Turbulent components of surface currents
are relatively low around this structure (not shown,
see http://pelabox.univ-lr.fr:3838/pelagis/DdSeaByc/;
Tab: '1. Data'), except at its most northern part. Again,
the location of this typical frontal structure con-
curred with the location of hot spots of fitted MI as
well as with the location of bycaught common dol-
phin strandings, mainly occurring on the coast of the
south Finistère region during the summer mortality
season (http://seamap.env.duke.edu).

4.1.3. Other drivers of seasonal differences

The distribution of common dolphins in the Bay
of Biscay differs between seasons (Laran et al. 2017).
In winter, they are more numerous on the shelf,
whereas in summer they are mainly distributed on
the shelf break. This probably results in spatial dif-
ferences in the local abundance of common dolphin

between seasons, directly affecting bycatch risk.
Seasonal differences in bycatch-induced mortality
may further be linked to annual biological cycles of
common dolphins and fish. The calving and mating
season of short-beaked common dolphins of the
northeast Atlantic is from April to September, with a
likely peak of activity in July and August (Murphy et
al. 2005, 2009). Feeding strategies and movement
patterns of dolphins differ during those periods, re-
sulting in a different response to oceanographic pro-
cesses. This could also be true with respect to dol-
phins' preys, many of which present an annual
biological cycle. As for the fisheries, they might like-
wise change their searching strategies, or their fish-
ing gear, depending on the season and target spe-
cies. The fishing effort and fishing targeted areas
also likely differ between seasons (ICES 2020). A
combination of those parameters might affect the dif-
ference between the 2 seasons regarding patterns of
oceanographic variables effects on bycatch mortality.

Although the mechanisms were different between
winter and summer mortality seasons, an influence of
oceanographic processes and especially seasonal
thermal front structures was highlighted for the 2
seasons of mortality. The association between marine
predator and thermal fronts was shown and studied
elsewhere (Doniol-Valcroze et al. 2007, Gannier &
Praca 2007, Cox et al. 2018, Scales et al. 2018), and
thermal fronts were identified as priority conserva-
tion areas for marine megafauna (Scales et al. 2014).
If correlations are complex and not straightforward,
this study still provides first hints of an association
between short-beaked common dolphins and this
type of mesoscale frontal structure in the Bay of Bis-
cay, likely linked to bycatch mortality.

4.2. Between-month variations

For all 3 oceanographic covariates, in winter, inter-
pretation is further complicated as their effects on the
mortality index often switched sign from one month
to the next. From a statistical point of view, there
could be 2 reasons for this. First, because we con-
ducted different models for each month, there is no
dependence between consecutive months. This could
be improved by putting additional model structure
such as a first order random walk for example. How-
ever, this would also substantially increase the com-
putational burden, multiplying the number of obser-
vations (3487 spatial cells 'observed' once per month
and per year, so 334 752 observations for the com-
plete dataset). Second, it is likely that the effect of
oceanographic processes on the mortality of dolphins due to bycatch is not linear or additive. Non-linear relationships with environmental covariates are routinely evidenced in marine species distribution modelling (Forney 2000, Stoner et al. 2001, Tew Kai & Marsac 2010) and fisheries (Maury et al. 2001, Walsh & Kleiber 2001, Zagaglia et al. 2004) or even bycatch studies (Žydelis et al. 2011, Di Tullio et al. 2015, Hahlbeck et al. 2017). Given the extra complexity of our response variable, it is quite conceivable that the influence of oceanographic features on short-beaked common dolphins bycatch is not linear. Combined with the fact that oceanographic processes might have different effects depending on their location and timing of appearance (that might vary from month to month), it could explain changes in signs of mean monthly coefficients of the oceanographic covariates. Again, relaxing the linear assumption would be possible with, for example, functional regression, but is beyond the scope of the present study.

There could also be ecological reasons for these between month variations. The winter oceanographic dynamics of the Bay of Biscay is characterized by high-frequency processes with rapid changes, sometimes within a month. Changes in the sign of OVs effects on the MI could be due to dynamic processes affecting common dolphins and fisheries interactions. The possibility that the effect of OVs could change from one month to the next motivated the use of random slopes. Taking thermal fronts as an example: these can aggregate preys of both common dolphins and predatory fish and do not necessarily increase bycatch risk. If fronts are quite wide and persistent, they could be wide enough so that they could attract fisheries on, say, their shore side and common dolphin on their seaside: spatial overlap between common dolphins and fisheries would be limited. However, if thermal fronts are very spatially localized, then spatial overlap between common dolphin and fisheries could be substantial, and bycatch risk would be higher. Yet, the frontal dynamics in the Bay of Biscay shows seasonal and interannual variability (Yelekçi et al. 2017). The width and intensity of thermal fronts can therefore vary from one month to the next (and for a given month, from one year to the next), which may also explain fickle statistical relationships.

### 4.3. The bycatch mortality of 2019: a predictable event?

Bycatch mortality in 2019 was the highest ever observed on the French Atlantic coast (Peltier et al.

2019), exceeding 2017 levels that were also unprecedented at that time. Still, crossvalidation results suggested it was statistically predictable as the observed monthly mortality tally from the dataset was within the confidence interval of the prediction (Fig. 7). Furthermore, results of the repetition scenarios showed how 2019's mortality was possibly due to circumstances already met in the past, either related to oceanographic conditions, or to other component included in models through the spatial process. The MI in 2019 was more accurately predicted when 2019 was considered as a repetition of 2013 and 2017 (Table 1, Fig. 7). For RS 2019 = 2017, the observed monthly mortality of 2019 was close to predicted values from oceanographic covariates (in 2019) and past relationships with OVs, both for winter and summer mortality seasons. As mentioned above, bycatch mortality in 2017 was also unprecedented at that time as an all-high record. 2017 indeed stands out as a singular year regarding the effect of all 3 covariates on bycatch risk (Fig. 5). This suggests maybe singular oceanographic phenomena in 2017 that were especially conductive of bycatch events, phenomena which recurred in 2019. A more surprising result was the remarkable good prediction made for winter months, when 2019 was considered as a repetition of 2013 (RS 2019 = 2013). Observed mortality in 2013 from the dataset was 156 individuals (as fresher stranded carcasses showing bycatch evidence or stranded during an intense mortality event, see Section 2), whereas it was 317 in 2017 and 413 in 2019, through January, February and March. However, the observed mortality in January 2013 was close to that of January 2017 (Fig. 3) and was the second highest value for this month in the study period. In addition, mean coefficient estimates of *eke* and *sst* were part of the highest from 2012 to 2018. In January, the mean coefficient estimate of *mean_sstGrad* in 2013 was close to the mean coefficient estimate in 2017. In January, the mean effect of the mean temperature gradient was also the same in 2013 and in 2017. Despite the lower levels of observed mortality in 2013, it seems that winter oceanographic conditions in 2013, 2017 and 2019 favored common dolphin and fishery interactions, with an intensifying trend over the years. Given the complexity of our results regarding oceanographic covariate effects on bycatch risk (with inter-annual and between-month variations), the additional value of repetition scenarios was to highlight whether similar associations occurred in different years.

Overall, these results showed that, from a statistical point of view, with few oceanographic covariates,
<figure><img src="image_165.png" /><figcaption>Fig. 7. Prediction for 2019. Monthly total mean mortality index from the dataset ('Observed MI') and from fitted values ('Predicted MI') for 2019 resulting from the 12 spatiotemporal models inferred with integrated nested Laplace approximations (INLA), with (1) a cross-validation (CV) procedure, i.e. the mortality index was predicted from previous years data and from values of oceanographic covariates for 2019 (2) repetition scenarios (RS), analogous to the latter except 2019 was considered as a repetition of a previous year of the study period. Grey dotted lines: 95%CI of fitted values</figcaption></figure>

the unprecedented bycatch mortality of 2019 was at the very least unsurprising as the pattern observed in 2019 could be accurately predicted from 2012–2018 data. This highlights how gaining knowledge about environmental influences on interactions between short-beaked common dolphins and fisheries could have great conservation and management value (Scales et al. 2018). In this study, we adopted a modelling approach to determine whether bycatch mortality can be explained by oceanographic processes. However, with evidenced oceanographic influences on bycatch incidence, a similar analysis could be designed to forecast potential hotspots of mortality. The modelling approach would have to be different, as explanatory and predictive modelling are distinct and are associated with different ‘modelling paths’ (Shmueli 2010). For example, other sources of environmental data, available in a smaller timeframe, might have to be considered for near real-time prediction. Remote sensing could be an option, as it is already used in the support of fisheries activities (Santos 2000). The forecasting of high mortality zones might make it possible to define dynamic time-area closures for fisheries (Cosandey-Godin et al. 2015), which represent a promising tool to reduce bycatch pressure on the common dolphin population while limiting the economic impact (Dunn et al. 2011, Maxwell et al. 2015, Hazen et al. 2018). In this regard, further studies modelling the dynamic distribution of

areas of co-occurrence of both short-beaked common dolphins and fisheries in the Bay of Biscay would be highly valuable for the conservation of this population.

## 4.4. Limitations and prospects for improvement

Several limitations can be imposed regarding the present work and results. First of all, strandings give only minimal estimations of bycatch-caused mortality. Some carcasses may never strand, and some others strand in advanced states of decomposition that do not allow their examination nor the identification of the cause of death of the animals. Moreover, due to the reverse drift modelling, the MI is associated with an uncertainty that is difficult to measure completely. Peltier et al. (2012) quantified the average distance between observed stranding locations of animals that were tagged when released from fishing boats and stranding locations predicted with MOTHY (27.1 ±

<table><caption>Table 1. Root mean squared error (RMSE) for cross-validation (CV) and repetition scenarios (RS)</caption><thead><tr><th></th><th>CV</th><th colspan="8">RS</th></tr><tr><th></th><th>2019</th><th>2012</th><th>2013</th><th>2014</th><th>2015</th><th>2016</th><th>2017</th><th>2018</th></tr></thead><tbody><tr><th>RMSE</th><td>48.4</td><td>51.2</td><td>32.8</td><td>43.6</td><td>54.2</td><td>45.0</td><td>30.1</td><td>48.1</td></tr></tbody></table>
24.5 km), which is accurate enough for practical purposes. Yet, there is no estimate of uncertainty associated with drift points as no information is available on where dolphins died. Developing the reverse modelling approach, on which the MI depends, was precisely motivated by inferring high-probability at-sea mortality areas (Peltier & Ridoux 2015, Peltier et al. 2020).

The spatial scale for the present analysis (approximately 89 km⁻²) was coarse compared to that of OVs. This coarse resolution required upscaling for OVs, which could have smoothed out too many high frequency processes such as fronts and eddies. It is thus possible that the relationship between the MI and OVs was weakened, impacting the results concerning oceanographic coefficients and explanatory power. The association between the MI and oceanographic variables is likely scale dependent as are associations between oceanographic covariates and marine top predators' distribution and foraging (Logerwell & Hargreaves 1996, Fauchald et al. 2000, Pinaud & Weimerskirch 2005, Pirotta et al. 2014). The link of interest here could be detectable at a specific spatial scale and not at another, and for high frequency processes the scale could be relatively small. Getting mortality data at a finer scale (and defining the spatial uncertainty associated with the reverse drift modelling) would then improve our understanding of the physical and ecological processes at stake.

Lastly, the models we developed for this study accounted for a low proportion of variance of the MI (25% at the highest, in February). This limits the range for interpretation, especially for months when this proportion was below 10%. This was expected: we aimed at quantifying links that are, at best, indirect between OVs and a marine top predator. Many trophic levels might be involved in the relationships of interest here, from phytoplankton to zooplankton and fish, each of these levels being influenced by several environmental and physiological factors. Moreover, the response variable we used is only a proxy of realized bycatch risk. Nonetheless, models performed well enough to reproduce annual patterns of total monthly MI from 2012 to 2018 (Appendix). Models performed better for the winter season with increased mortality and to a lesser extent for the summer mortality season, which were of particular interest regarding the seasonal pattern of bycatch. The CV results finally confirmed the acceptable performance of our models, suggesting an association between chosen oceanographic variables and bycatch-induced mortality of common dolphins in the Bay of Biscay.

This association should not be over-interpreted: relationships highlighted here are not straightforward, and they have been highlighted only for a few months, and not for all years. The limited explanatory power of the oceanographic covariates even for months when bycatch mortality was relatively high suggests that bycatch risk in the Bay of Biscay may be influenced by other non-environmental factors, e.g. characteristics of fisheries like the type of fishing gear used, the overall fishing effort or the timeline of the fishing operation (day or night). There were no significant changes in the pattern of the relationship between oceanographic processes and bycatch risk that could help explain the rise in bycatch observed from 2016 onward. This further suggests an influence of additional factors not considered in the present study. This rise of bycatch could be associated with a change in the distribution of common dolphins in the Bay of Biscay or a change in fishing practices. Bycatch risk is likely affected by a mix of different factors, oceanographic processes being one among others.

## 5. CONCLUSIONS

This study demonstrated the relevance of standing data and modelled OVs from coastal operational oceanography to better understand this intense bycatch phenomenon which causes a serious threat to the short-beaked common dolphin population of the northeast Atlantic (Peltier et al. 2016, Murphy et al. 2019). The need to explicitly integrate spatial and temporal dimensions for such complex interactions was justified: our results effectively highlighted the spatiotemporal influences of structuring oceanographic processes on the risk of common dolphin bycatch in the Bay of Biscay. However, the underlying mechanisms remain unclear, and the relationship between environmental characteristics and dolphin bycatch must not be overinterpreted, as our model only indirectly links OVs with dolphin bycatch. Because of the dynamic aspects of both structuring oceanographic processes and dolphin and fishery distributions, further research focusing on smaller time scales is needed. For instance, extreme standing events are observed on a weekly timeframe (Peltier et al. 2020), suggesting short periods of intense interactions between dolphins and fisheries, so a similar analysis on a weekly timeframe and focusing on these events might be more informative. Furthermore, we showed that the pattern of seasonal bycatch events can be directly inferred from past data on strandings and oceanographic conditions. This
highlights the value of establishing functional relationships to eventually identify adapted and acceptable management measures to ensure both the long-term conservation of the population of short-beaked common dolphins of the Bay of Biscay and fisheries activities.

Acknowledgements. We warmly thank all members of the French and UK stranding networks for their continuous effort in collecting stranding data. We thank Marie-Pierre Etienne and Etienne Rivot for their insightful comments on our work. We also thank the community participating in the development of INLA, and especially Havard Rue for his devotion towards helping INLA users. Finally, we thank the reviewers for their insightful comments which helped to significantly improve the quality of the manuscript. The Observatoire PELAGIS is funded by the ministry in charge of the environment, and by Communauté d'Agglomération de la Ville de La Rochelle, with support from Région Poitou-Charentes and the European Union. EOVs were computed by the Shom and partially funded by the ministry in charge of the environment under the Marine Strategy Framework Directive (MSFD) monitoring program.

LITERATURE CITED

<figure><img src="image_166.png" /></figure>

Bailey H, Thompson P (2010) Effect of oceanographic features on fine-scale foraging movements of bottlenose dolphins. Mar Ecol Prog Ser 418:223-233

Bakun A (1996) Patterns in the ocean: ocean processes and marine population dynamics. California Sea Grant, La Jolla, CA

<figure><img src="image_167.png" /></figure>

Bakun A (2006) Fronts and eddies as key structures in the habitat of marine fish larvae: opportunity, adaptive response and competitive advantage. Sci Mar 70(S2):105-122

<figure><img src="image_168.png" /></figure>

Barlow DR, Bernard KS, Escobar-Flores P, Palacios DM, Torres LG (2020) Links in the trophic chain: modeling functional relationships between *in situ* oceanography, krill, and blue whale distribution under different oceangoing regimes. Mar Ecol Prog Ser 642:207-225

<figure><img src="image_169.png" /></figure>

Becker EA, Forney KA, Ferguson MC, Foley DG, Smith RC, Barlow J, Redfern JV (2010) Comparing California current cetacean-habitat models developed using *in situ* and remotely sensed sea surface temperature data. Mar Ecol Prog Ser 413:163-183

Besag J, Kooperberg C (1995) On conditional and intrinsic autoregressions. Biometrika 82:733-746

<figure><img src="image_170.png" /></figure>

Besag J, York J, Mollié A (1991) A Bayesian image restoration with two applications in spatial statistics. Ann Inst Stat Math 43:1-59

<figure><img src="image_171.png" /></figure>

Best BD, Halpin PN, Read AJ, Fujioka E and others (2012) Online cetacean habitat modeling system for the US East Coast and Gulf of Mexico. Endang Species Res 18:1-15

<figure><img src="image_172.png" /></figure>

Blangiardo M, Cameletti M, Baio G, Rue H (2013) Spatial and spatio-temporal models with R-INLA. Spat Spatiotemporal Epidemiol 4:33-49

<figure><img src="image_173.png" /></figure>

Bleck R (2002) An oceanic general circulation model framed in hybrid isopycnic-Cartesian coordinates. Ocean Model 4:55-88

<figure><img src="image_174.png" /></figure>

Brownell RL Jr, Reeves RR, Read AJ, Smith BD and others (2019) Bycatch in gillnet fisheries threatens Critically

Endangered small cetaceans and other aquatic megafauna. Endang Species Res 40:285-296

<figure><img src="image_175.png" /></figure>

Castro J, Couto A, Borges FO, Cid A, Laborde MI, Pearson HC, Rosa R (2020) Oceanographic determinants of the abundance of common dolphins (*Delphinus delphis*) in the south of Portugal. Oceans 1:165-173

<figure><img src="image_176.png" /></figure>

Cosandey-Godin A, Krainski ET, Worm B, Flemming JM (2015) Applying Bayesian spatiotemporal models to fisheries bycatch in the Canadian Arctic. Can J Fish Aquat Sci 72:186-197

<figure><img src="image_177.png" /></figure>

Cox SL, Embling CB, Hosegood PJ, Votier SC, Ingram SN (2018) Oceanographic drivers of marine mammal and seabird habitat-use across shelf-seas: a guide to key features and recommendations for future research and conservation management. Estuar Coast Shelf Sci 212: 294-310

<figure><img src="image_178.png" /></figure>

Croll DA, Tershy BR, Hewitt RP, Demer DA and others (1998) An integrated approach to the foraging ecology of marine birds and mammals. Deep Sea Res II 45:1353-1371

<figure><img src="image_179.png" /></figure>

D'Agrosa C, Lennert-Cody CE, Vidal O (2000) Vaquita bycatch in Mexico's artisanal gillnet fisheries: driving a small population to extinction. Conserv Biol 14:1110-1119

<figure><img src="image_180.png" /></figure>

Daniel P, Jan G, Cabioc'h F, Landau Y, Loiseau E (2002) Drift modeling of cargo containers. Spill Sci Technol Bull 7: 279-288

Dars C, Dabin W, Demaret F, Meheust E and others (2019) Les échouages de mammifères marins sur le littoral français en 2019. Rapport scientifique de l'Observatoire Pelagis, La Rochelle Université et CNRS

<figure><img src="image_181.png" /></figure>

Di Tullio JC, Fruet PF, Secchi ER (2015) Identifying critical areas to reduce bycatch of coastal common bottlenose dolphins *Tursiops truncatus* in artisanal fisheries of the subtropical western South Atlantic. Endang Species Res 29:35-50

<figure><img src="image_182.png" /></figure>

Diaz López B, Methion S, Giralt Paradell O (2019) Living on the edge: overlap between a marine predator's habitat use and fisheries in the northeast Atlantic waters (NW Spain). Prog Oceanogr 175:115-123

<figure><img src="image_183.png" /></figure>

Doniol-Valcroze T, Berteaux D, Larouche P, Sears R (2007) Influence of thermal fronts on habitat selection by four rorqual whale species in the Gulf of St. Lawrence. Mar Ecol Prog Ser 335:207-216

<figure><img src="image_184.png" /></figure>

Dunn DC, Boustany AM, Halpin PN (2011) Spatio-temporal management of fisheries to reduce by-catch and increase fishing selectivity. Fish Fish 12:110-119

Fauchald P, Erikstad KE, Skarsfjord H (2000) Scale-dependent predator-prey interactions: the hierarchical spatial distribution of seabirds and prey. Ecology 81:773-783

Fioravanti G, Martino S, Cameletti M, Cattani G (2021) Spatio-temporal modelling of PM10 daily concentrations in Italy using the SPDE approach. Atmos Environ 248:118192

<figure><img src="image_185.png" /></figure>

Forney KA (2000) Environmental models of cetacean abundance: reducing uncertainty in population trends. Conserv Biol 14:1271-1286

<figure><img src="image_186.png" /></figure>

Gabry J, Simpson D, Vehtari A, Betancourt M, Gelman A (2019) Visualization in Bayesian workflow. J R Stat Soc Ser A Stat Soc 182:389-402

<figure><img src="image_187.png" /></figure>

Gannier A, Praca E (2007) SST fronts and the summer sperm whale distribution in the northwest Mediterranean Sea. J Mar Biol Assoc UK 87:187-193

<figure><img src="image_188.png" /></figure>

García-Barón I, Santos MB, Saavedra C, Astarloa A, Valeiras J, Garcia Barcelona S, Louzao M (2020) Essential ocean variables and high value biodiversity areas: targets for the conservation of marine megafauna. Ecol Indic 117:106504
<figure><img src="image_189.png" /></figure>

Gardner B, Sullivan PJ, Epperly S, Morreale SJ (2008) Hierarchical modeling of bycatch rates of sea turtles in the western North Atlantic. Endang Species Res 5:279-289

<figure><img src="image_190.png" /></figure>

Gelman A, Hwang J, Vehtari A (2014) Understanding predictive information criteria for Bayesian models. Stat Comput 24:997-1016

<figure><img src="image_191.png" /></figure>

Hahlbeck N, Scales KL, Dewar H, Maxwell SM, Bograd SJ, Hazen EL (2017) Oceanographic determinants of ocean sunfish (Mola mola) and bluefin tuna (Thunnus orientalis) bycatch patterns in the California large mesh drift gillnet fishery. Fish Res 191:154-163

<figure><img src="image_192.png" /></figure>

Hazen EL, Scales KL, Maxwell SM, Briscoe DK and others (2018) A dynamic ocean management tool to reduce bycatch and support sustainable fisheries. Sci Adv 4:eaar3001

<figure><img src="image_193.png" /></figure>

Hughes KM, Dransfeld L, Johnson MP (2014) Changes in the spatial distribution of spawning activity by north-east Atlantic mackerel in warming seas: 1977-2010. Mar Biol 161:2563-2576

<figure><img src="image_194.png" /></figure>

Huret M, Bourriau P, Doray M, Gohin F, Petitgas P (2018) Survey timing vs. ecosystem scheduling: degree-days to underpin observed interannual variability in marine ecosystems. Prog Oceanogr 166:30-40

<figure><img src="image_195.png" /></figure>

Hyrenbach DK, Forney KA, Dayton PK (2000) Viewpoint: marine protected areas and ocean basin management. Aquat Conserv 10:437-458

ICES (2020) Workshop on fisheries emergency measures to minimize bycatch of short-beaked common dolphins in the Bay of Biscay and harbour porpoise in the Baltic Sea (WKEMBYC). ICES Sci Rep 2:354

<figure><img src="image_196.png" /></figure>

Ijsseldijk LL, ten Doeschate MT, Brownlow A, Davison NJ and others (2020) Spatiotemporal mortality and demographic trends in a small cetacean: strandings to inform conservation management. Biol Conserv 249:108733

<figure><img src="image_197.png" /></figure>

Jaramillo-Legorreta A, Rojas-Bracho L, Brownell RL Jr, Read AJ, Reeves RR, Ralls K, Taylor BL (2007) Saving the vaquita: immediate action, not more data. Conserv Biol 21:1653-1655

<figure><img src="image_198.png" /></figure>

Jaramillo-Legorreta AM, Cardenas-Hinojosa G, Nieto-Garcia E, Rojas-Bracho L and others (2019) Decline towards extinction of Mexico's vaquita porpoise (Phocoena sinus). R Soc Open Sci 6:190598

<figure><img src="image_199.png" /></figure>

Johnson PC (2014) Extension of Nakagawa & Schielzeth's R2GLMM to random slopes models. Methods Ecol Evol 5:944-946

Kuiken T (1994) Diagnosis of bycatch in cetaceans. Proceedings of the Second ECS Workshop on Cetacean Pathology. European Cetacean Society, Montpellier

<figure><img src="image_200.png" /></figure>

Lahaye S, Gouillon F, Baraille R, Pichon A, Pineau-Guillou L, Morel Y (2011) A numerical scheme for modeling tidal wetting and drying. J Geophys Res Oceans 116:1-20

<figure><img src="image_201.png" /></figure>

Lambert C, Virgili A, Pettex E, Delavenne J, Toison V, Blanck A, Ridoux V (2017) Habitat modelling predictions highlight seasonal relevance of marine protected areas for marine megafauna. Deep Sea Res II 141:262-274

Lanz E, López-Martínez J, Nevárez-Martínez M, Dworak JA (2009) Small pelagic fish catches in the Gulf of California associated with sea surface temperature and chlorophyll. CCOFI Rep 50:134-146

<figure><img src="image_202.png" /></figure>

Laran S, Authier M, Blanck A, Doremus G and others (2017) Seasonal distribution and abundance of cetaceans within French waters—Part II. The Bay of Biscay and the English Channel. Deep Sea Res II 141:31-40

<figure><img src="image_203.png" /></figure>

Legendre P, Dale MR, Fortin MJ, Gurevitch J, Hohn M, Myers D (2002) The consequences of spatial structure for

the design and analysis of ecological field surveys. Ecol-ogy 25:601-615

<figure><img src="image_204.png" /></figure>

Levin SA (1992) The problem of pattern and scale in ecol-ogy. Ecology 73:1943-1967

<figure><img src="image_205.png" /></figure>

Lewison RL, Crowder LB, Read AJ, Freeman SA (2004) Understanding impacts of fisheries bycatch on marine megafauna. Trends Ecol Evol 19:598-604

<figure><img src="image_206.png" /></figure>

Logerwell EA, Hargreaves NB (1996) The distribution of sea birds relative to their fish prey off Vancouver Island: opposing results at large and small spatial scales. Fish Oceanogr 5:163-175

<figure><img src="image_207.png" /></figure>

Mannocci L, Dabin W, Augeraud-Véron E, Dupuy JF, Barbraud C, Ridoux V (2012) Assessing the impact of bycatch on dolphin populations: the case of the common dolphin in the Eastern North Atlantic. PLOS ONE 7: e32615

<figure><img src="image_208.png" /></figure>

Martinez-Minaya J, Cameletti M, Conesa D, Pennino MG (2018) Species distribution modeling: a statistical review with focus in spatio-temporal issues. Stochastic Environ Res Risk Assess 32:3227-3244

<figure><img src="image_209.png" /></figure>

Maury O, Gascuel D, Marsac F, Fonteneau A, Rosa ALD (2001) Hierarchical interpretation of nonlinear relationships linking yellowfin tuna (Thunnus albacares) distribution to the environment in the Atlantic Ocean. Can J Fish Aquat Sci 58:458-469

<figure><img src="image_210.png" /></figure>

Maxwell SM, Hazen EL, Lewison RL, Dunn DC and others (2015) Dynamic ocean management: defining and conceptualizing real-time management of the ocean. Mar Policy 58:42-50

<figure><img src="image_211.png" /></figure>

Mellin C, Russell BD, Connell SD, Brook BW, Fordham DA (2012) Geographic range determinants of two commercially important marine molluscs. Divers Distrib 18: 133-146

<figure><img src="image_212.png" /></figure>

Miller PI, Scales KL, Ingram SN, Southall EJ, Sims DW (2015) Basking sharks and oceanographic fronts: quantifying associations in the north-east Atlantic. Funct Ecol 29:1099-1109

<figure><img src="image_213.png" /></figure>

Moore MJ, Mitchell GH, Rowles TK, Early G (2020) Dead cetacean? Beach, bloat, float, sink. Front Mar Sci 7:333

Moraga P (2019) Geospatial health data: modeling and visualization with R-INLA and Shiny. CRC Press, Boca Raton, FL

<figure><img src="image_214.png" /></figure>

Morel Y, Baraille R, Pichon A (2008) Time splitting and linear stability of the slow part of the barotropic component. Ocean Model 23:73-81

<figure><img src="image_215.png" /></figure>

Murphy S, Collet A, Rogan E (2005) Mating strategy in the male common dolphin (Delphinus delphis): what gonadal analysis tells us. J Mammal 86:1247-1258

<figure><img src="image_216.png" /></figure>

Murphy S, Winship A, Dabin W, Jepson PD and others (2009) Importance of biological parameters in assessing the status of Delphinus delphis. Mar Ecol Prog Ser 388: 273-291

Murphy S, Evans PG, Pinn E, Pierce GJ (2019) Conservation management of common dolphins: lessons learned from the north-east Atlantic. Aquat Conserv 131:137-166

<figure><img src="image_217.png" /></figure>

Murray KT, Orphanides CD (2013) Estimating the risk of loggerhead turtle Caretta caretta bycatch in the US mid-Atlantic using fishery-independent and -dependent data. Mar Ecol Prog Ser 477:259-270

<figure><img src="image_218.png" /></figure>

Nel DC, Lutjeharms JRE, Pakhomov EA, Ansorge IJ, Ryan PG, Klages NTW (2001) Exploitation of mesoscale oceanographic features by grey-headed albatross Thalassarche chrysostoma in the southern Indian Ocean. Mar Ecol Prog Ser 217:15-26
* Norris KS, McFarland WN (1958) A new harbor porpoise of the genus *Phocoena* from the Gulf of California. *J Mammal* 39:22

* Pebesma E (2018) Simple features for R: standardized support for spatial vector data. *R J* 10:439–446

* Peltier H, Ridoux V (2015) Marine megavertebrates adrift: a framework for the interpretation of stranding data in perspective of the European Marine Strategy Framework Directive and other regional agreements. *Environ Sci Policy* 54:240–247

* Peltier H, Dabin W, Daniel P, Van Canneyt O, Dorémus G, Huon M, Ridoux V (2012) The significance of stranding data as indicators of cetacean populations at sea: modelling the drift of cetacean carcasses. *Ecol Indic* 18:278–290

* Peltier H, Baagøe HJ, Camphuysen KC, Czeck R and others (2013) The stranding anomaly as population indicator: the case of harbour porpoise *Phocoena phocoena* in north-western Europe. *PLOS ONE* 8:e62180

* Peltier H, Jepson PD, Dabin W, Deaille R, Daniel P, Van Canneyt O, Ridoux V (2014) The contribution of stranding data to monitoring and conservation strategies for cetaceans: developing spatially explicit mortality indicators for common dolphins (*Delphinus delphis*) in the eastern North-Atlantic. *Ecol Indic* 39:203–214

* Peltier H, Authier M, Deaille R, Dabin W and others (2016) Small cetacean bycatch as estimated from stranding schemes: the common dolphin case in the northeast Atlantic. *Environ Sci Policy* 63:7–18

Peltier H, Authier M, Caurant F, Dabin W, Dars C and others (2019) Etat des connaissances sur les captures accidentelles de dauphins communs dans le golfe de Gascogne – Synthèse 2019. Rapport scientifique dans le cadre de la convention avec le MTES. Technical report, Observatoire PELAGIS—UMS 3462, La Rochelle Université/CNRS

* Peltier H, Authier M, Dabin W, Dars C and others (2020) Can modelling the drift of bycaught dolphin stranded carcasses help identify involved fisheries? An exploratory study. *Glob Ecol Conserv* 21:e00843

* Peltier H, Authier M, Caurant F, Dabin W and others (2021) In the wrong place at the wrong time: identifying spatiotemporal co-occurrence of bycaught common dolphins and fisheries in the Bay of Biscay (NE Atlantic) from 2010 to 2019. *Front Mar Sci* 8:617342

* Petitgas P, Doray M, Huret M, Masse J, Wouillez M (2014) Modelling the variability in fish spatial distributions over time with empirical orthogonal functions: anchovy in the Bay of Biscay. *ICES J Mar Sci* 71:2379–2389

* Pichon A, Correard S (2006) Internal tides modelling in the Bay of Biscay. Comparisons with observations. *Sci Mar* 70(S1):65–88

* Pichon A, Morel Y, Baraille R, Quaresma LS (2013) Internal tide interactions in the Bay of Biscay: observations and modelling. *J Mar Syst* 109–110(Suppl):S26–S44

* Pinaud D, Weimerskirch H (2005) Scale-dependent habitat use in a long-ranging central place predator. *J Anim Ecol* 74:852–863

* Pirotta E, Thompson PM, Miller PI, Brookes KL and others (2014) Scale-dependent foraging ecology of a marine top predator modelled using passive acoustic data. *Funct Ecol* 28:206–217

* Planque B, Bellier E, Lazoné P (2007) Modelling potential spawning habitat of sardine (*Sardina pilchardus*) and anchovy (*Engraulis encrasicolus*) in the Bay of Biscay. *Fish Oceanogr* 16:16–30

* Politikos DV, Huret M, Petitgas P (2015) A coupled movement and bioenergetics model to explore the spawning migration of anchovy in the Bay of Biscay. *Ecol Model* 313:212–222

R Core Team (2020) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna

Reeves RR, Berggren P, Crespo EA, Gales N and others (2005) Global priorities for reduction of cetacean bycatch. Technical report. World Wildlife Fund

* Rue H, Martino S, Chopin N (2009) Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. *J R Stat Soc Series B Stat Methodol* 71:319–392

* Santos AMP (2000) Fisheries oceanography using satellite and airborne remote sensing methods: a review. *Fish Res* 49:1–20

* Scales KL, Miller PI, Hawkes LA, Ingram SN, Sims DW, Votier SC (2014) On the front line: frontal zones as priority at-sea conservation areas for mobile marine vertebrates. *J Appl Ecol* 51:1575–1583

* Scales KL, Hazen EL, Jacox MG, Castruccio F, Maxwell SM, Lewison RL, Bograd SJ (2018) Fisheries bycatch risk to marine megafauna is intensified in Lagrangian coherent structures. *Proc Natl Acad Sci USA* 115:7362–7367

* Shmueli G (2010) To explain or to predict? *Stat Sci* 25:289–310

* Snyder S, Franks PJS, Talley LD, Xu Y, Kohin S (2017) Crossing the line: tunas actively exploit submesoscale fronts to enhance foraging success. *Limnol Oceanogr Lett* 2:187–194

* Spitz J, Chouvelon T, Cardinaud M, Kostecki C, Lorance P (2013) Prey preferences of adult sea bass *Dicentrarchus labrax* in the northeastern Atlantic: implications for bycatch of common dolphin *Delphinus delphis*. *ICES J Mar Sci* 70:452–461

* Stephenson F, Goetz K, Sharp BR, Mouton TL and others (2020) Modelling the spatial distribution of cetaceans in New Zealand waters. *Divers Distrib* 26:495–516

* Stoner AW, Manderson JP, Pessutti JP (2001) Spatially explicit analysis of estuarine habitat for juvenile winter flounder: combining generalized additive models and geographic information systems. *Mar Ecol Prog Ser* 213:253–271

* Tew Kai E, Marsac F (2010) Influence of mesoscale eddies on spatial structuring of top predators' communities in the Mozambique Channel. *Prog Oceanogr* 86:214–223

* Tew-Kai E, Quilfen V, Cachera M, Boutet M (2020) Dynamic coastal-shelf seascapes to support marine policies using operational coastal oceanography: the French example. *J Mar Sci Eng* 8:585

* Tobler WR (1970) A computer movie simulating urban growth in the Detroit region. *Econ Geogr* 46:234

* Turvey ST, Pitman RL, Taylor BL, Barlow J and others (2007) First human-caused extinction of a cetacean species? *Biol Lett* 3:537–540

* Valcu M, Kempenaers B (2010) Spatial autocorrelation: an overlooked concept in behavioral ecology. *Behav Ecol* 21:902–905

* Walsh WA, Kleiber P (2001) Generalized additive model and regression tree analyses of blue shark (*Prionace glauca*) catch rates by the Hawaii-based commercial longline fishery. *Fish Res* 53:115–131

Watanabe S (2010) Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. *J Mach Learn Res* 11:3571–3594

* Weltz K, Kock AA, Winker H, Attwood C, Sikweyiya M (2013) The influence of environmental variables on the
presence of white sharks, *Carcharodon carcharias* at two
popular Cape Town bathing beaches: a generalized
additive mixed model. PLOS ONE 8:e68554

Yelekçi Ö (2017) Submesoscale dynamics in the Bay of Biscay continental. PhD, Université Pierre et Marie Curie – Paris VI

Yelekçi Ö, Charria G, Capet X, Reverdin G, Sudre J, Yahia H (2017) Spatial and seasonal distributions of frontal activity over the French continental shelf in the Bay of Biscay. Cont Shelf Res 144:65–79

Zagaglia CR, Lorenzzetti JA, Stech JL (2004) Remote sensing data and longline catches of yellowfin tuna (*Thunnus albacares*) in the equatorial Atlantic. Remote Sens Environment 93:267–281

Zuur A, Ieno EN, Smith GM (2007) Analyzing ecological data, Chaps 6 & 7. Springer, New York, NY

Žydelis R, Lewison RL, Shaffer SA, Moore JE and others (2011) Dynamic habitat models: using telemetry data to project fisheries bycatch. Proc R Soci B Biol Sci 278: 3191–3200

Appendix.

<figure><img src="image_219.png" /></figure>

Fig. A1. Temporal series of the monthly total mean mortality index from the dataset ('observed', in blue) and from fitted values ('fitted', in purple) resulting from the 12 spatiotemporal models inferred with integrated nested Laplace approximations (INLA), from 2012 to 2018. Grey dotted lines: sum of the 0.025 (lower bound) and 0.975 (upper bound) quantiles of fitted values
<figure><img src="image_220.png" /></figure>

Check for
updates

# Estimating Cetacean Bycatch From Non-representative Samples (I): A Simulation Study With Regularized Multilevel Regression and Post-stratification

Matthieu Authier¹,²*, Etienne Rouby³ and Kelly Macleod⁴

¹Observatoire Pelagis, UMS 3462, CNRS-La Rochelle Université, La Rochelle, France, ²ADERA, Pessac, France, ³Centre d'Études Biologiques de Chizé, UMR 7372 CNRS-La Rochelle Université, Villiers-en-Bois, France, ⁴Joint Nature Conservation Committee, Aberdeen, United Kingdom

## OPEN ACCESS

**Edited by:**
Jeff Moore
Southwest Fisheries Science Center
(NOAA), United States

**Reviewed by:**
Peter Dillingham,
University of Otago, New Zealand
Devin S. Johnson,
Pacific Islands Fisheries Science
Center, National Marine Fisheries
Service, United States

**\*Correspondence:**
Matthieu Authier
mauthier@univ-lr.fr

**Specialty section:**
This article was submitted to
Marine Megafauna,
a section of the journal
Frontiers in Marine Science

**Received:** 03 June 2021
**Accepted:** 17 September 2021
**Published:** 20 October 2021

**Citation:**
Authier M, Rouby E and Macleod K (2021) Estimating Cetacean Bycatch From Non-representative Samples (I): A Simulation Study With Regularized Multilevel Regression and Post-stratification. *Front. Mar. Sci.* 8:719956.
doi: 10.3389/fmars.2021.719956

Bycatch, the non-intentional capture or killing of non-target species in commercial or recreational fisheries, is a world wide threat to protected, endangered or threatened species (PETS) of marine megafauna. Obtaining accurate bycatch estimates of PETS is challenging: the only data available may come from non-dedicated schemes, and may not be representative of the whole fisheries effort. We investigated, with simulated data, a model-based approach for estimating PETS bycatch from non-representative samples. We leveraged recent development in the statistical analysis of surveys, namely regularized multilevel regression with post-stratification, to infer total bycatch under realistic scenarios of data sampling such as under-sampling or over-sampling when PETS bycatch risk is high. Post-stratification is a survey technique to re-align the sample with the population and addresses the problem of non-representative samples. Post-stratification requires to sub-divide a population of interest into potentially hundreds of cells corresponding to the cross-classification of important attributes. Multilevel regression accommodate this data structure, and the statistical technique of regularization can be used to predict for each of these hundreds of cells. We illustrated these statistical ideas by modeling bycatch risk for each week within a year with as few as a handful of observed PETS bycatch events. The model-based approach led to improvements, under mild assumptions, both in terms of accuracy and precision of estimates and was more robust to non-representative samples compared to more design-based methods currently in use. In our simulations, there was no detrimental effects of using the model-based even when sampling was representative. Estimating PETS bycatch ideally requires dedicated observer schemes and adequate coverage of fisheries effort. We showed how a model-based approach combining sparse data typical of PETS bycatch and recent methodological developments can help when both dedicated observer schemes and adequate coverage are challenging to implement.

**Keywords:** marine megafauna, conservation, modeling, statistics, Bayesian, anthropogenic removals
# 1. INTRODUCTION

Bycatch, the non-intentional capture or killing of non-target species in commercial or recreational fisheries, is a world wide threat to protected, endangered or threatened species (PETS) of marine megafauna (Gray and Kennelly, 2018), including seabirds (Martin et al., 2019), elasmobranchs (Pacoureau et al., 2021) and cetaceans (Avila et al., 2018). Bycatch in fishing gears, such as gillnets, is currently driving some small cetacean species to extinction (Brownell et al., 2019; Jaramillo-Legorreta et al., 2019). The European Commission recently issued infringement procedures against several Member States for failing to correctly transpose some provisions of European environmental law (the Habitats Directive, Council Directive 92/43/EEC), in particular the obligations related to the establishment of a coherent monitoring scheme of cetacean bycatch<fnref n="1" />. The Data Collection Framework (DCF) provides a common framework in the European Union (EU) to collect, manage, and share data within the fisheries sector (Anonymous, 2019a). The Framework indicates that the Commission shall establish a Multi-Annual Union Programme (EU-MAP) for the collection and management of fisheries data which should be inclusive of data that allows the assessment of fisheries' impact on marine ecosystems. With respect to PETS (including cetaceans), the collection of high quality data usually requires a dedicated sampling scheme and methodology, and is generally different from those applied under the DCF (Stransky and Sala, 2019): "EU MAP remains not well-suited for the dedicated monitoring of rare and protected bycatch in high-risk fisheries since its main focus is the statistically-sound random sampling of all commercial fisheries (Ulrich and Doerner, 2021, p. 126)." In practice, the introduction of any programme on PETS bycatch under the DCF may be met with caution because of its perceived potential to disrupt data collection for fisheries management (Stransky and Sala, 2019). This perception implicitly relegates PETS bycatch as a side issue for fishery management rather than an integral part of it. It may explain the usually poor quality of bycatch data on PETS (ICES, 2020a).

Recent EU legislation (Regulation 2019/1241), referred to as the Technical Measures Regulation (TMR), requires Member States to collect scientific data on cetacean bycatch for the following métiers: pelagic trawls (single and pair), bottom-set gillnets and entangling nets; and high-opening trawls (Anonymous, 2019b). Unlike its predecessor (Council Regulation EC No. 812/2004), this Regulation does not require the establishment of dedicated observer schemes for cetacean bycatch data collection (Dolman et al., 2020). Furthermore, only vessels of an overall length of 15 m or more are to be monitored, but these represent a small fraction of the European fleet (less than 10% in 2019)<fnref n="2" />. This vessel length criterion introduces bias in the bycatch monitoring data as the sample of vessels larger than 15 m is almost certainly dissimilar to smaller

vessels. Even within the sample of vessels that are monitored, pragmatic considerations can complicate sampling. For example, in the United States, observer sampling trips are allocated first by region, port, and month, then randomly to vessels of particular categories within those monthly and spatial strata (ICES, 2009). Random allocation of observers to vessels follows sound statistical methodology and increases the likelihood of collecting unbiased data (Babcock and Pikitch, 2003). In France, observer days are allocated by port and by month for each fishery, but the exact vessel allocation is then negotiated and left at the discretion of skippers (ICES, 2009). Allocation is no longer random as skippers may only accept observers when cetacean bycatch risk is low (Benoît and Allard, 2009). Non-random allocation means potential bias in the collected data for monitoring bycatch as the sub-sample of skippers accepting an observer may be very different from skippers refusing to do so (Babcock and Pikitch, 2003).

One pragmatic solution bypassing observers is to mandate skippers to self-declare the non-intentional capture or killing of any PETS, as already required under the DCF (Anonymous, 2019a). In France, a national law from 2011 mandate fisheries to declare (without fear of prosecution) the bycatch of any cetacean species, but this law remained largely unknown to French fishermen until late 2019 (Cloâtre, 2020). In general, self-reported PETS bycatch data are sub-optimal as they may be heavily biased, non-representative (ICES, 2009) and typically provide poor information on which to base management decisions (National Marine Fisheries Service, 2004). Once again, the set of skippers who choose to declare bycatch may differ markedly from those who do not: for example the former take the extra time required to fill logbooks and thus provide accurate data while the latter do not. If this behavior is correlated to other attributes, e.g., a more acute awareness of threats to PETS resulting in practices that tend to minimize impact on PETS, data collected from skippers reporting bycatch would not be representative. There may also be an element of skippers genuinely forgetting to log PETS bycatch in the bustle of the fishing operation but this is random and unlikely to introduce bias. In addition, ground-truthing, for example with remote-electronic monitoring (REM; Course et al., 2020), would be required in order to ensure the quality and accuracy of self-reported data before their statistical analyses.

Another hurdle, of the statistical kind, with cetacean bycatch is the low frequency of these events. Assuming that implementing a representative sampling program were feasible, if bycatch is a rare event (Komoroske and Lewison, 2015), then few events would be observed for realistic sampling effort (Babcock and Pikitch, 2003; ICES, 2009). This paucity of observed event means a large uncertainty in statistical estimates: with a bycatch rate of the order of 0.01 event per fishing operation, a sample size of 1,100 observed operations would be required to obtain, in the best case scenario (no bias, statistical independence, etc.), the US recommended coefficient of variation of 30% (National Marine Fisheries Service, 2005, 2016; ICES, 2009; Carretta and Moore, 2014). The amount of observer coverage needed to reach this precision depends on fishery size and trip duration (Babcock and Pikitch, 2003). In practice, the sampling error depends on the overall design of the survey, of which the sample size

<footnote n="1">¹https://ec.europa.eu/info/news/july-infringements-package-commission-moves-against-member-states-not-respecting-eu-energy-rules-2019-jul-26_en</footnote>

<footnote n="2">²https://appsso.eurostat.ec.europa.eu/nui/show.do?dataset=fish_fleet_alts&lang=en</footnote>
is only one factor: for example a larger sample size could be
needed if there are large “skipper-effects” as the same vessels
would contribute fishing operations, and these would not be
statistically independent. With a small sample size, uncertainty
may be so large as to prevent using estimates altogether, even if
one were to assume no bias in the data (Babcock and Pikitch,
2003). Given this challenge and the lack of uptake of dedicated
monitoring programmes of cetacean bycatch in Europe over the
last decade or more (Sala et al., 2019), it would appear prudent
to seek methods of analysis that can handle the few and non-
representative data available to robustly estimate bycatch rates.

The problem of having non-representative samples to carry
out statistical analyses is ancient (Hansen and Hurwitz, 1946)
and widespread: it pops up in many applied disciplines, including
election forecasting (Wang et al., 2015; Kiewiet De Jonge et al.,
2018), political sciences (Lax and Phillips, 2009; Zahorski, 2020),
social sciences (Halsny, 2020), addiction studies (Rhem et al.,
2020) or epidemiology (Zhang et al., 2014; Downes et al., 2018).
In these disciplines, there are also intrinsic limits on improving
the representativeness of sampling. For example, in polling, non-
response rates can be above 90% (Forsberg, 2020). In other
cases, some populations of interest may be hard to reach (Rhem
et al., 2020), or answers may not be honest (St. John et al.,
2014). Challenges lie in the accurate estimation of quantities
of scientific interest (e.g., the true magnitude of bycatch in
a fishery; Babcock and Pikitch, 2003) with the construction
of statistical weights that can calibrate a non-representative
survey sample to the population targets. Such weights are
implicit with simple random sampling where each unit in a
population has the same, non-nil, probability of being included
in the sample. When inclusion probabilities differ between units,
weights inversely proportional to the former can be used to adjust
the sample. However, constructing survey weights is in general
more elaborate than using inverse probabilities of selection in the
sample (Gelman, 2007). Model-based approaches, and multilevel
regression modeling with post-stratification in particular, has
become an attractive alternative to weighting to adjust non-
representative samples (Gelman, 2007).

Multilevel regression modeling allows researchers to
summarize how predictions of an outcome of scientific interest
vary across statistical units defined by a set of attributes or
covariates (Gelman et al., 2021, p. 4): for example bycatch events
are a binary outcome at the fishing operation level (a unit)
associated with attributes, such as date-time, location, gears and
vessels (e.g., Palka and Rossman, 2001). Post-stratification is a
standard technique to generalize inferences from a sample to the
population by adjusting for known discrepancies between the
former and the latter. Post-stratification is a form of adjustment
whereby statistical units are sorted out according to an auxiliary
variable (hereafter a stratum) after completion of data collection;
stratum-level effects (i.e., effects within each stratum or cell) are
then estimated, and finally averaged with weights proportional
to stratum size to obtain the population-level estimate. Post-
stratification differs from blocking as the latter is done before
data collection to ensure balance and representativeness at the
design stage. Post-stratification is a *post hoc* statistical adjustment
done at the analysis stage: it can remove bias, but at the price

of an increased variance of estimates. Lennert et al. (1994)
provided an early example of model-based estimates of bycatch
with post-stratification.

In small samples post-stratification can degrade estimate
precision, especially if the number of strata is large as each
stratum will typically include very few data, or even not a single
datum (the so-called “small-area” problem). In practice, adequate
post-stratification may require handling hundreds of cells (the
crossing of several attributes; e.g., week by statistical area by
gears). Some predictions for each cell may be too noisy, especially
if there are sparse or no data for that particular combination
of attributes. Multilevel regression can offer a solution as it
borrows strength from similar units to improve and stabilize
(i.e., regularize) these predictions (Cam, 2012). In other words,
multilevel regression allows an efficient use of a sparse sample
to estimate the outcome of interest within each cell, even if
these cells are very numerous (e.g., several hundreds). The key
insight of combining multilevel regression modeling with post-
stratification is thus: even if observations are not a representative
sample of the population of interest, it may be possible to
construct a regression model to first predict unobserved cases,
and then post-stratify to average the fitted regression model’s
predictions over the population of interest (Gelman et al., 2021,
p. 313). Good predictions may be obtained with regularization
by means of multilevel models with structured priors (including
so called “random-effects” models). The latter can increase
precision by inducing shrinkage of parameter estimates across
similar post-stratification cells, where similarity is encoded in the
model specification (e.g., by using random effects that assume
exchangeability). The amount of shrinkage, or partial-pooling
across cells, is model-based and thus data-driven. However,
in order to be able to leverage the information in the data,
some model structure on the parameters of interest is necessary
hence the need for structured priors. Relying on a model rather
than just empirical means of the response variable addresses
the bias-variance problem intrinsic to having a large number
of cells in post-stratification, and leverages the large toolbox of
regression-based models.

Technically, when data arise as signal plus noise, overfitting
occurs when a regression model captures too much of the
noise compared to the signal; that is in using an ill-conditioned
(unstable) model that will provide an excellent in-sample fit
but make poor out-of-sample predictions (Authier et al., 2017b;
George and Ročková, 2021). Overfitting may result when using
richly parametrized models without using adequate estimation
methods such as regularization to stabilize parameter estimates
and buffer them against noise (Gelman et al., 2021, p. 459–460).
Weakly-informative priors in a Bayesian framework regularize
the estimation of the large number of parameters that may be
present in a multilevel model. Multilevel modeling takes into
account complex data structures with structured prior models
for batches of parameters; the simplest example are so-called
“random effects” whereby a common (Gaussian) distribution
centered on zero and with an unknown variance to be estimated
for data is assumed for a group of parameters; for example
years or sites (Cam, 2012). This common distribution for the
parameters is a prior model, and this model for parameters means
that the latter are not independently estimated but in tandem
according to the postulated prior model. For example, Sims et al.
(2008) used a model-based approach to obtain spatially smoothed
estimates of bycatch in a gillnet fishery. Spatial-smoothing (also
known as "small-area estimation"; Fay and Herriot, 1979) was
used to stabilize estimated bycatch rates by using a Conditional
Autoregressive prior model that leverages information from
spatial neighbors to improve the prediction at a specific location.
Prior models add some soft constraints to the overall model
and these constraints are very useful in data sparse settings to
mitigate variance and bias in predictions. In other words, these
prior models represent additional assumptions about the data,
assumptions, which if approximately correct, add information in
the analyses and increase the precision and stability of predictions
at the cost of a usually small estimation bias. Introducing bias
to reduce variance is a common statistical technique known as
shrinkage or regularization (George and Ročková, 2021).

Regularized multilevel regression with post-stratification is thus the combination of several important ideas to obtain accurate predictions (Gao et al., 2019). First, post-stratification is a survey technique to re-align the sample with the population and addresses the problem of non-representative samples. In practice, post-stratification requires to sub-divide the population of interest into many cells corresponding to the combination of important attributes. Multilevel regression can be used to accommodate all these cells in a single model, but the problem has now moved to how to obtain useful estimates for all these cells, which can number in the several hundreds. Regularization solves this estimation problem: it introduces model-driven bias in statistical estimates in order to stabilize them. These new developments in the statistical analysis of non-representative samples may help in obtaining a better quantification of bycatch rates and numbers. Our aim is to assess with simulations, the potential of regularized multilevel regression with post-stratification for analyzing already collected bycatch data, with the full knowledge that these data are non-representative and biased in several respects. These biases in sampling are manifold (see above): bias may be due to regulation exempting certain vessels (e.g., no monitoring for vessels smaller than 15 m); to non-dedicated observers or because sampling is driven for other purposes than bycatch monitoring of PETS (commercial discards, stock assessment); or in the case of dedicated schemes, to over-sampling a few "cooperative" skippers or focusing sampling in métiers with the highest or lowest bycatch risk. Our focus will be narrower, honing in on specific sampling scenarios whereby observer coverage is correlated to bycatch risk. In other words, we will assess the potential of regularized multilevel regression with post-stratification to estimate accurately bycatch numbers with samples preferentially collected either during low- or high-bycatch risk periods. Our investigation is largely framed from our knowledge on small cetacean bycatch in European waters, such as short-beaked common dolphin (*Delphinus delphis*, lower observer coverage when bycatch risk is higher) in the Bay of Biscay (Peltier et al., 2021) or harbor porpoises (*Phocoena phocoena*, higher observer coverage when bycatch risk is higher) in the Celtic Seas (Tregenza et al., 1997). In the remainder, we first introduce methods and notations to

detail the proposed model to perform multilevel regression
with post-stratification with bycatch data, using dolphins as
an example. Next, we explain our data simulation scenarios
and how we emulate non-representative sampling. We then
compare the results (i.e., estimates of bycatch) from the proposed
modeling approach with those from the method currently used
by the working group on bycatch of protected species from
the International Council for the Exploration of the Sea (ICES
WGBYC) before concluding on some recommendations for
future investigations.

## 2. MATERIALS AND METHODS

We carried out Monte Carlo simulations to assess the ability
of regularized multilevel regression with post-stratification to
estimate bycatch risk and bycatch numbers from representative
and non-representative samples. ICES WGBYC collate data
through an annual call from dedicated and DCF surveys
collecting data on the bycatch of PETS through onboard
observers or REM. These surveys may be qualified as "design-
based" in the sense that, ideally, a representative coverage of
fisheries would be sought in order to scale up the observed sample
to the whole population using ratio-estimators. There are many
caveats around the use of these ratio-estimators as EU MAP is not
well-suited for monitoring PETS bycatch (Ulrich and Doerner,
2021). Given these shortcomings in the collection of bycatch data
under EU MAP, the data available to ICES WGBYC are unlikely
to be representative of fisheries of interest but nevertheless, ratio-
estimators are used as part of a Bycatch Risk Approach (BRA)
to identify relative risk of bycatch across species and metiers
(ICES, 2018). Cetacean bycatch observer programmes may aim
at achieving a pre-specified precision for bycatch rates (with a
coefficient of variations less than 30%; National Marine Fisheries
Service, 2005, 2016; ICES, 2009; Carretta and Moore, 2014).
Achieving this is very difficult in practice, and a given coverage
of effort deployed by the total fleet is, instead, aimed at: for
example 10% (5%) for pair-trawlers (level-3 métier PTM) larger
(smaller) than 15 m in France. Data from onboard observer
programmes are then used to estimate total bycatch using ratio
estimators (Lennert et al., 1994; Julian and Beeson, 1998; Amandé
et al., 2012) and the bootstrap or a classical approach (Clopper-
Pearson) for uncertainty quantification (ICES, 2018, p. 57). We
used an approach similar to that of WGBYC (hereafter referred
to as a "design-based" approach) as a benchmark to compare
against results from regularized multilevel regression with post-
stratification. We honed in on the accurate estimation of the
number of bycatch events for a complete fleet. We assume that
information on the total effort deployed by a fleet operating
in a spatial domain are available and measured without error.
This assumption is necessary to scale estimates from the sample
to the population. We also assumed that there are no false-
negatives in the sample, that is no bycatch event went unrecorded
by onboard observers (assuming thereby a dedicated observer
programme). These two assumptions are customary with ratio
estimators, whether design- or model-based, and do not deviate
from current norms. We assume however that these population
data on total effort can be disaggregated at a finer temporal scale in order to post-stratify on calendar weeks. This assumption of accurate measurement of effort at the week-level is crucial for post-stratification.

## 2.1. Notations

The logit transform maps a quantity $p \in ]0, 1[$ to the real line: $\logit(p) = \log(\frac{p}{1-p})$. Its inverse is denoted by $\logit^{-1}(x) = \frac{1}{1+e^{-x}}$ (sometimes called the “expit” transform). Let $y_{ijkl}$ denote the ith fishing operation of vessel j in week k of year l, with $y_{ijkl} = 1$ if a bycatch event occurs and 0 otherwise:

$$y_{ijkl} \sim \text{Bernoulli} \left( p_{ijkl} = \logit^{-1} \left( \mu + \beta_{kl} + \alpha_j \right) \right) \quad (1)$$

where $p_{ijkl}$ is the product of the probability of a bycatch event occurring and the probability of dolphin presence. This unconditional probability $p_{ijkl}$, or “bycatch risk” hereafter, is not indexed by i: although there may be several fishing operations of vessel j in week k of year l, the risk is assumed constant over these. Bycatch risk is a function of several parameters (on a logit scale): $\mu$ is the intercept (overall risk), $\alpha_j \sim N(0, \sigma_{vessel})$ are (unstructured, normal random effects) vessel-effects accounting for heterogeneity (e.g., “fishing style” of skippers); and $\beta_{kl}$ are time effects, modeled with a Gaussian Process. A Gaussian process is written as $GP(m, c)$ where m and c are the mean and covariance functions respectively (Gelman et al., 2014, p. 501). The Gaussian Process prior on the vector of week effects in year l, $\beta_l$, defines this vector as a random function for which the values at any week 1,...,k,...,w are drawn from a w-multivariate normal distribution:

$$\beta_{1l}, \dots, \beta_{wl} \sim N((m(1), \dots, m(w)), \Omega) \quad (2)$$

with mean m and covariance $\Omega$. The function c specifies the covariance between any 2 weeks k and k', with $\Omega$ an $w \times w$ covariance matrix with element $\Omega(k, k') = c(k, k')$. A Matérn covariance function of order $\frac{3}{2}$ and range parameter fixed to $\frac{3}{2}$ was assumed: $c(k, k') = \sigma_{year}^2 \times (1 + \frac{2\sqrt{3}\times d(k-k')}{3}) \times \exp(-\frac{2\sqrt{3}\times d(k-k')}{2})$, where $d(k-k')$ is the temporal distance (in weeks) between weeks k and k'. The distance function was the absolute difference between calendar weeks within the same year: $d(k-k') = |k-k'|$. The choice of the Matérn covariance function translate an assumption of smoothness in the temporal profile of bycatch risk: bycatch risk is assumed to change gradually across weeks, with no abrupt increase or decrease. The range parameter is fixed and not estimated from data. This choice represents an additional assumption whereby the temporal correlation is 0.05 after 4 weeks corresponding to temporal independence after a month. This choice is to some extent arbitrary and represents an additional assumption. In theory, the range parameter could also be estimated from data but we assumed a data sparse setting with limited information (more so with Bernoulli data) to estimate this parameter.

The mean function m of the Gaussian process was modeled (on a logit scale) with a first order random walk, which was

evaluated at specific values $k \in [1,...,w]$ corresponding to week number within a year:

$$\left\{ \begin{array}{ll} (m(1), \dots, m(w)) = (\varepsilon_1, \dots, \varepsilon_w) & k = 1 \\ \varepsilon_k = 0 & k > 1 \\ \varepsilon_{k+1} \sim N(\varepsilon_k, \sigma_{\text{week}}) & k > 1 \end{array} \right. \quad (3)$$

The order of the random walk prior was assumed fixed at 1 and not estimated from data. This prior choice smooths the first order differences between adjacent elements of $\varepsilon$ and represents an additional assumption, mainly to limit the number of parameters to estimate from the typically sparse data on bycatch. A random walk was chosen as an effective way to reveal the shape of the average risk profile without specifying a family of parametric curves.

The model in Equation (1) is a decomposition of bycatch risk into a time-varying component (at the week-scale, Equation 3; and with an interaction with year, Equation 2) and time-invariant component which can be interpreted as fishing-style effects whereby some skippers may have consistent practices that increase or decrease bycatch risk. Importantly, bycatch risk is modeled here with no attempt to model dolphin presence directly as relevant data to do so may be missing in the general case. Bycatch risk is thus to be estimated for each week of a year, and each of these weeks represent de facto a stratum. In any applied case, additional factors, such as statistical area, may need to be included in Equation (1) for improved realism. For simplicity, we did not consider space in simulations, and solely focused on time.

## 2.2. Data Simulation

To test the ability of model 1 to estimate bycatch risk, data were simulated (Figure 1).

1. Bycatch probability conditional on dolphin presence was constant and set to 0.3, that is roughly one fishing operation out of 3 generates a bycatch event when dolphins are present (corresponding to a high risk fishery, e.g., the trawl fishery in the Bay of Biscay).

2. Dolphin presence is seasonal (loosely inspired from the observed pattern of common dolphin in the Bay of Biscay where abundance is higher closer to the coasts in winter; Laran et al., 2017): it peaks at the beginning and end of the year, but quickly drops to 0 for roughly 2 thirds of a year.

3. A fishery of 20 vessels is operating all year round, with an overall activity rate of 80% each week (that is, for any week, $20 \times \frac{80}{100} = 16$ vessels are fishing). Each fishing day (5 days per week), on average 2.3 fishing operations are carried out. The expected total number of fishing operations for a year is $5 \times 52 \times 2.3 \times 16 \approx 10,000$. These values were loosely taken from an exploratory analysis of onboard observer data collected on PTM flying the French flag. During each of these operations, a bycatch event may occur depending on dolphin presence at the time and on a skipper-specific risk factor (drawn randomly from a normal distribution with scale parameter set to $\frac{\log(2)}{3}$ to induce moderate heterogeneity on a logit scale; Authier et al., 2017a).
<figure><img src="image_221.png" /></figure>

4. Observers are accepted onboard vessels either with a constant probability of 0.05 corresponding to a coverage of 5% of all fishing operations (unbiased sampling scenario) or with a probability that covaries with dolphin presence (biased sampling scenarios). In the latter case, realized coverage is a random variable. With under-sampling, the bulk of the observer data is collected when bycatch risk (the product of dolphin presence and bycatch probability) is nil (Figure 1). With over-sampling, the bulk of the observer data is collected when bycatch risk is high but no data are collected when the risk is nil (Figure 1).

5. In a year, the number of fishing operations is ≈ 10,000, and the number of bycatch events ≈ 300, which yields a rate of ≈ 3%. This rate is not large, but is not extremely rare either.

Bycatch events were simulated for each fishing operations during a day when an observer was present from a Bernoulli distribution according to the product of bycatch probability given dolphin presence and dolphin presence probability for that day. If no observer was present, no data were recorded. The data-generating mechanism used a parametric function for dolphin presence probability and was different from the statistical model used to analyze the data (see https://gitlab.univ-lr.fr/mauthier/regularized_bycatch). For each sampling scenario, 100 datasets were generated for 1, 5, 10, or 15 years. All data simulations were carried out in R v4.0.1 (R Core Team, 2020). When simulating only 1 year of data, Equation (2) is not necessary as there is no between-year variation to estimate: the model can be simplified with the omission of $β_l$. Our Monte Carlo study had a comprehensive factorial design crossing (a) sampling regime
(either unbiased or not) and (b) sample size as controlled with the number of years for which the observer programme was assumed to have been in operation.

## 2.3. Estimation

Estimation of the parameters of model 1 from simulated data was carried out in a Bayesian framework using programming language Stan (Carpenter et al., 2017) called from R v.4.0.1 (R Core Team, 2020) with library Rstan (Stan Development Team, 2020). Stan uses Hamiltonian dynamics in Markov chain Monte Carlo (MCMC) to sample values from the joint posterior distribution (Carpenter et al., 2017). Weakly-informative priors

$$ \text{were used for regularization: } \begin{cases} \mu \sim \mathcal{N}(0, \frac{3}{2}) \\ \text{prop} \sim \mathcal{D}(1, 1, 1) \\ \sigma_{\text{total}} \sim \mathcal{GG}(\frac{1}{2}, \frac{1}{2}, \frac{\log 2}{10}) \end{cases} $$

where $\mathcal{D}$() denotes the Dirichlet distribution for modeling proportions (such that $\sum_{i=1}^3 \text{prop}_i = 1$) and $\mathcal{GG}$() the Gamma-Gamma distribution for scale parameters (Griffin and Brown, 2017; Pérez et al., 2017). With this simplex parametrization, chosen to improve mixing and ease estimation with Monte Carlo methods (He et al., 2007), the several variance components of

$$ \text{the model were: } \begin{cases} \sigma_{\text{vessel}}^2 = \sigma_{\text{total}}^2 \times \text{prop}_1 \\ \sigma_{\text{week}}^2 = \sigma_{\text{total}}^2 \times \text{prop}_2 \\ \sigma_{\text{year}}^2 = \sigma_{\text{total}}^2 \times \text{prop}_3 \end{cases} $$

These priors are weakly-informative (Gabry et al., 2019); the prior for the intercept covers the whole interval between 0 and 1 on the probability scale but is informative on the logit scale. The prior for the scale (square-root of the variance) is heavy tailed and has a median set to $\frac{\log 10}{2}$ (Griffin and Brown, 2017; Pérez et al., 2017), which translate an assumption about the plausible range of variations in bycatch risk spanning a priori two full order of magnitude from one tenth to a ten-fold increase compared to the mean bycatch rate. Thirty random realizations from our choice of priors are depicted on Figure 2: the whole interval between 0 and 1 is covered, and between-week variations can be large or small.

For each simulated dataset, four chains were initialized from diffuse random starting points (Carpenter et al., 2017, p. 20) and run for a total of 1,000 iterations, discarding the first 500 as warm-up. Default settings for the No-U-Turn Sampler (NUTS) were changed to 0.99 for adapt delta and 15 for max treedepth (Hoffman and Gelman, 2014). NUTS uses Hamiltonian Dynamics in MCMC and typically requires shorter runs than other MCMC algorithms both to reach convergence and to obtain an equivalent Effective Sample Size from the posterior (Hoffman and Gelman, 2014; Monnahan et al., 2017). Parameter convergence was assessed using the $\hat{R}$ statistics (Vehtari et al., 2019) and assumed if $\hat{R} < 1.025$. Upon diagnosing convergence of all parameters, a combined sample of 4 × 500 = 2,000 MCMC values were obtained to approximate the joint posterior distribution. Let $\hat{\mu}^{[m]}$, $\hat{\beta}_{kl}^{[m]}$, $\hat{\sigma}_{\text{vessel}}^{[m]}$ denote the mth MCMC sample for parameters $\mu$, $\beta_{kl}$ and $\sigma_{\text{vessel}}$. Bycatch risk $\hat{P}_{j^*kl}$ for a randomly chosen vessel $j^*$ operating in week k of year l was computed from the mth MCMC draw from the joint posterior

distribution as:

$$ \hat{P}_{j^*kl}^{[m]} = \logit^{-1} \left( \hat{\mu}^{[m]} + \hat{\beta}_{kl}^{[m]} + \hat{\alpha}_{j^*}^{[m]} \right) \quad (4) $$

where $\hat{\alpha}_{j^*}^{[m]} \sim \mathcal{N}(0, \hat{\sigma}_{\text{vessel}}^{[m]})$. This predicted bycatch risk incorporates between-vessel variability, that is it takes into account the fishing style of skippers. The predicted risk (on a logit scale) for a random chosen skipper is $\hat{\alpha}_{j^*}^{[m]}$ and was drawn from the posterior predictive distribution: not all skippers may be observed in the sample, and but the subset of skippers that accept an observer can be used to estimate a between-skipper variance in bycatch risk. In practice, the number of fishing operations carried out in the course of a week in a year by individual skippers is unknown, although the aggregated number of fishing operations may be known. If totals by skippers were available, and all skippers had been sampled, it would be more efficient to use skipper-specific estimated risk, but we did not assume that this would necessarily be the case.

The total number of bycatch events, $T_{\text{bycatch}}$ was estimated as the average over the 2,000 MCMC draws from the posterior:

$$ \hat{T}_{\text{model-based bycatch}} = \frac{1}{2000} \sum_{m=1}^{2000} \left( \sum_{l=1}^{n_{\text{year}}} \sum_{k=1}^{n_{\text{week}}} \hat{P}_{j^*kl}^{[m]} \times N_{kl} \right) \quad (5) $$

where $N_{kl}$ is the total number of fishing operations that took place is week k of year l. The total number of strata for post-stratification was $n_{\text{year}} \times n_{\text{week}}$, with a maximum of 15 × 52 = 780 cells. Highest Posterior Density credible intervals at the 80% level were computed with function HPDinterval from package coda (Plummer et al., 2006) for uncertainty evaluation. Equation (5) is an instance of a ratio-estimator with post-stratification, except that it uses model-based estimates of bycatch risk. This model-based approach regularizes estimates with partial pooling (Gelman and Shalizi, 2013): the variance of estimates is greatly reduced by introducing some bias with structured priors (Gao et al., 2019). Our results were benchmarked against an approach similar to that of ICES WGBYC whereby total number of bycatch events was estimated as:

$$ \hat{T}_{\text{design-based bycatch}} = \sum_{l=1}^{n_{\text{year}}} \left( \tilde{p}_l \times \sum_{k=1}^{n_{\text{week}}} N_{kl} \right) \quad (6) $$

where $\tilde{p}_l$ is the average bycatch risk estimated as the mean from the observed sample in year l. Confidence intervals at the 95% level were computed using either the bootstrap or the Clopper-Pearson approach as customary in ICES WGBYC. Both were considered as the Clopper-Pearson approach is known for being more conservative: it produces confidence intervals that above the nominal level (i.e., wider than necessary) but generates non-nil confidence intervals even if no bycatch has been observed (Northridge et al., 2019). In practice, ICES WGBYC often pooled several years to stabilize the estimate of $\tilde{p}$ (e.g., ICES 2018, p. 57–58; Carretta and Moore, 2014): Equation (6) translate an ideal case that is rarely met in practice.
<figure><img src="image_222.png" /></figure>

ICES WGBYC usually works on bycatch rates (in number of PETS per unit effort), not bycatch risk. We focused on risk for simplicity, but scaling bycatch risk to a rate is straightforward by multiplying with the average number of PETS bycaught in a bycatch event. Dolphin presence was seasonal in the data-generating mechanism for simulations: pitching a method that can explicitly accommodate such seasonality against one that does not may be viewed as knocking down a strawman. However, current estimates of PETS bycatch in Europe are stratified by flag. ICES statistical areas, and métiers but not by season (e.g., Table 2 p. 17 in ICES 2019; Northridge et al., 2019, p. 27). The comparison remains relevant and topical as it matches current practices.

# 3. RESULTS

Convergence across all simulations and scenarios was assumed to be reached, with all $\hat{R} < 1.025$, for all parameters. For each simulation, chains were combined in a single sample of 2,000

values to approximate the joint posterior distribution of the model defined by Equations (1), (2), and (3).

## 3.1. Design- vs. Model-Based Approach

Comparing the design- and model-based approach was done with simulating 1 year of data. When data sampling was unbiased, both the design- and model-based approach were able to recover the true number of bycatch events (Figure 3; Table 1). Estimates of bycatch events were statistically unbiased but their precision low with a (frequentist 95%) confidence or (Bayesian 80%) credible interval (CI) as large as 100% of the point estimate (Table 1), as could be expected with only 15 bycatch events were recorded on average by onboard observers (Table 1). With under-sampling, design-based estimates were negatively biased (that is, they were under-estimates) whereas model-based estimates were still unbiased on average (Figure 3; Table 1). With over-sampling, design-based estimates were positively biased (that is, they were over-estimates) but so were model-based estimates, although bias was 5 times smaller (Figure 3; Table 1). In all cases, coverage was 100% but largely as a result of low precision:
<figure><img src="image_223.png" /><figcaption>TABLE 1 | Statistical properties of estimates from the design- and model-based approach.</figcaption></figure>

<table><thead><tr><th>Method</th><th>Uncertainty</th><th>Data sampling</th><th>n<sub>years</sub></th><th>Bias (%)</th><th>Coverage (%)</th><th>Width of CI (%)</th><th>n<sub>obs</sub></th></tr></thead><tbody><tr><td>Design-based</td><td>Bootstrap</td><td>Unbiased</td><td>1</td><td>3.5</td><td>100.0</td><td>102.5</td><td>15</td></tr><tr><td>Design-based</td><td>Clopper-Pearson</td><td>Unbiased</td><td>1</td><td>3.5</td><td>100.0</td><td>115.0</td><td>15</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Unbiased</td><td>1</td><td>3.6</td><td>100.0</td><td>120.4</td><td>15</td></tr><tr><td>Design-based</td><td>Bootstrap</td><td>Under-</td><td>1</td><td>-83.5</td><td>100.0</td><td>195.0</td><td>5</td></tr><tr><td>Design-based</td><td>Clopper-Pearson</td><td>Under-</td><td>1</td><td>-83.5</td><td>100.0</td><td>259.6</td><td>5</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Under-</td><td>1</td><td>3.0</td><td>100.0</td><td>204.3</td><td>5</td></tr><tr><td>Design-based</td><td>Bootstrap</td><td>Over-</td><td>1</td><td>121.0</td><td>100.0</td><td>46.1</td><td>63</td></tr><tr><td>Design-based</td><td>Clopper-Pearson</td><td>Over-</td><td>1</td><td>121.0</td><td>100.0</td><td>50.1</td><td>63</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Over-</td><td>1</td><td>22.1</td><td>100.0</td><td>78.6</td><td>63</td></tr></tbody></table>

One year of data was simulated 100 times. Bias of point estimate, coverage of frequentist 95% confidence or (Bayesian 80%) credible interval (CI) and precision (as CI width relative to the point estimate) are reported. The last column indicates the average number of bycatch events ($n_{obs} = E[\sum_{i=1}^{n} y_{ik}]$) that were recorded by onboard observers during data sampling.
precision was very low with CI spanning some 200% of the
point estimate for the unbiased and under-sampling scenarios.
This low precision was the result of having to work with as
few as 5 observed bycatch events on average (Table 1). Precision
improved with over-sampling, but was still as high as 50% of
the point (over-)estimate. The model-based approach was well-
calibrated in both the unbiased and under-sampling scenarios
(Figure 4): model-based estimates were on average equal to the
truth whereas this was only the case with design-based estimates
when sampling was unbiased. In addition, the model-based
approach was able to recover the temporal profile of bycatch risk
(Figure 5) in these two scenarios, but with an increased accuracy
and precision if sampling was unbiased. In the over-sampling
scenario, both the design- and model-based approaches were not
well -calibrated (Figure 4) and the model-based approach over-
estimated bycatch risk when no data were collected (Figures 1, 5).

## 3.2. Model-Based Approach With Several Years of Data

With several years of data, the model-based approach was able
to yield nearly unbiased estimates: the bias was smaller than
3 bycatch events when sampling was unbiased, but as large
as 10 (on average) with biased sampling and 3 years of data.
The precision of estimates improved with several years of data,
as expected with larger sample size. Precision of model-based
estimates with over-sampling were already acceptable with 3
years of data: an 80% credible interval width of 50% corresponds
to a coefficient of variation of $\frac{50}{2.5} \approx 20\%$ assuming a normal
distribution for the posterior. The model-based approach allowed
to obtain estimates at the weekly scale (Figure 6): these estimates
were approximately unbiased in the unbiased and over-sampling
scenarios, but were biased for the under-sampling scenario. In
that latter case, the bias was correlated with the temporal pattern
used to simulate dolphin presence (Figure 1): it was the largest
when dolphin presence was at its highest but positive at the
beginning of a year and negative at the end of the same year. Both
biases were greatly attenuated with increased sample size.

# 4. DISCUSSION

Using Monte-Carlo simulations, we investigated the statistical
properties of a model-based approach, regularized multilevel
regression with post-stratification, to estimate the total number
of bycatch events in a fishery operating year-round. Simulations
were broadly informed from the case of common dolphins and
pair-trawlers in the Bay of Biscay and from harbor porpoises
and set-gillnets in Celtic Seas. A salient feature of simulations
was biased sampling with observers being preferentially accepted
onboard when bycatch risk was either high or low. Data
simulations in that latter case, which is the most realistic one
in the Bay of Biscay (Peltier et al., 2016), resulted in as few
as 5 observed bycatch events per year on average (Tables 1, 2).
This aligns with the ubiquitous description of small cetacean
bycatch being a rarely observed event. It was nevertheless possible
to fit a regularized multilevel regression model on these data.
Importantly, estimates from this model-based approach were

statistically less biased than the design-based estimates when
sampling was biased. Model-based estimates were, however,
imprecise but this is largely to be expected (Amandè et al.,
2012), especially with as few as 5 observed bycatch events per
year. The design-based approach was also imprecise, even in
the unbiased data sampling scenario of 5% coverage of the
fleet, which is not reached in practice (Anonymous, 2016; ICES,
2020b). The design-based approach was very sensitive to how
data were collected: this approach severely under- or over-
estimated bycatch when sampling was biased, whereas the model-
based approach was still well-calibrated with under-sampling, but
not with over-sampling (Figure 4).

Biases in onboard observer data are pervasive and widely
acknowledged (Babcock and Pikitch, 2003; Benoît and Allard,
2009; Peltier et al., 2016). Enforcing coverage as required to
achieve a pre-specified precision in estimates can be challenging
in practice. For example, in 2016, France only achieved a
coverage rate less than 2% for most métiers and concluded on
the impossibility of scaling-up observed bycatch rates to the
whole fleet (Anonymous, 2016, p. 24). There were, however, 9
bycatch events of common dolphins in pair-trawlers targeting
European hake (Merluccius merluccius). From these numbers,
bycatch was described a “rare” event (Anonymous 2016, p.
23). Such a conclusion would be warranted if sampling were
representative, in which case the design-based estimate could be
used, even though its precision would still be very low. On the
other hand, with under-sampling, this conclusion is misleading
as our simulations further illustrated: although only 5 bycatch
events were observed on average (Figure 4), the true number
of bycatch events was on average 60 times larger (Figure 4).
In our simulations, the true bycatch rate was on average ≈ 3% over
a year, which is not rare, but not frequent either. Moreover,
interviews with French skippers deploying trawls or gillnets in
the Bay of Biscay revealed that more than 80% of respondents
declared to having experienced at least one small cetacean
bycatch event in a year (Cloâtre, 2020). Such a large proportion
contradicts the idea of common dolphin bycatch being a rare
event in the Bay of Biscay, but rather suggest severe biases
in onboard observer data that result in the rare reporting of
bycatch events, rather than a rarity of events *per se*. The common
dolphin in the Bay of Biscay illustrates how under-sampling
may distort the perception of bycatch as a very rare event
when it can, in fact, be widespread. This is a catch-22 situation
whereby cetacean bycatch is described as a rare event because
it is rarely reported, and this perceived rarity may serve to
argue against ambitious dedicated monitoring programmes out
of cost-effective considerations, thereby preventing to dispel the
initial misconception.

Finding an optimal sampling plan for fisheries with rare
bycatch events is a long standing problem (ICES, 2009). Several
strategies have been attempted: for example in the United
States, one strategy is “pulsed sampling” whereby a particular
fishery or métier is very heavily sampled for a short period
of time in order to maximize the chance for observers to
record any bycatch that might occur (ICES, 2009). This pulsed
sampling strategy corresponds to our over-sampling scenario
wherein monitoring effort is positively correlated with bycatch
<figure><img src="image_224.png" /></figure>

risk. Under this scenario, the absence of any sampling at all when bycatch risk was low was detrimental to the accurate estimation of bycatch events with our model. Model-based estimates were, however, less biased than design-based estimates. Arguably, this comparison is somewhat artificial as a correct comparison would use all the available information and uses estimators that are season-specific to account for under-sampling when bycatch risk is low if such a period is known to the investigator. Notwithstanding this shortcoming, model-based estimates represented an improvement and allowed to infer the bycatch risk profile accurately, especially with several years of data.

We showed with our Monte-Carlo simulations that regularized multilevel regression with post-stratification can nevertheless be used to analyze bycatch data despite concerns about non-representative sampling. Model-based approaches (Palka and Rossman, 2001), with post-stratification (Lennert

et al., 1994), or machine learning (Carretta et al., 2017), or multilevel regression (Sims et al., 2008; Martin et al., 2015) have previously been used to estimate bycatch rates. Traditional, design-based, ratio estimates are biased if sampling is biased; imprecise if observer coverage is low (as is the usual case in the North East Atlantic; see for example Figure 14, p. 114 in ICES, 2020b); and volatile if bycatch events are only observed occasionally (Carretta et al., 2017). The traditional remedy to stabilize estimates and improve precision is to bypass year-specific estimation and pool several years together (Carretta and Moore, 2014; ICES, 2018). This pragmatic solution improves precision but does not address the problem of biased sampling. It also introduces estimation bias for any year-specific estimates by pooling completely several years in order to stabilize the variance of estimates (ICES, 2009, p. 36): any between-year differences are thus ignored in order to obtain a better precision of estimates. It is a reasonable approach in practice, but one
<figure><img src="image_225.png" /></figure>

that can be improved. Model-based approaches offer a trade-off between no-pooling (keeping all years separate) and complete-pooling with a third option: partial pooling or regularization (Gelman and Shalizi, 2013). Regularization is a general term for statistical procedures that give more stable estimates. Our model-based approach achieves regularization by leveraging, via a structured prior model (Equations 2 and 3, see section 2), the within-year information at the weekly scale. The result were more stable and accurate annual bycatch estimates at the cost of some modeling assumptions and weakly-informative priors. Importantly, weekly estimates could also be obtained with our model-based approach.

Our model-based approach is semi-parametric as it uses a random walk prior to learn from the data the weekly pattern in bycatch risk. This prior is also ensuring some smoothness in the temporal risk profile as it translates an assumption on the correlation between 2 consecutive weeks. This random walk

model remains simple as the order is fixed to 1. We further expanded this model to allow for between-years variation in the weekly risk profile with a Gaussian Process prior (Neal, 1998; Goldin and Purse, 2016). Importantly, these two prior choices (a random walk and a Gaussian Process prior) add structure to the model and help in leveraging the information present in the sparse data typical of onboard observer programmes. Even when with over-sampling, these choices were not detrimental as model-based estimates were statistically unbiased and precise with 3 years of data (Table 2). The explicit consideration of time effects is key to mitigate bias in sampling. In our simulations, dolphin presence was caricaturally seasonal, and observers could be preferentially allowed on fishing vessels when dolphins were less or more likely to be present (Figure 1). Our model was still able to provide statistically unbiased estimates of bycatch in those scenarios, although these estimates were very imprecise with under-sampling. However, they were not
<figure><img src="image_226.png" /></figure>

more imprecise than the traditional (but biased) design-based estimates (Table 1) if 80% credible interval were used. In addition to being unbiased, these estimates could also reveal with accuracy the temporal risk profile (Figure 5). It is important to keep in mind here that our model is different from the data-generating model used in simulating data: our results were not simply an instance of using a true model, which is impossible in practice as a model is by definition a simplification used to capture the salient features of a phenomenon. Our model had some shortcomings: for example, bias increased with 3 years of data compared to 1 year for the under-sampling scenario (contrast Tables 1, 2). This increased bias (toward the prior model) was the result of partial pooling but came with a gain in precision as evidenced in the width of credible intervals. The bias progressively wore off with more years of data, illustrating thereby the attractiveness of partial pooling and structured priors to regularize estimates (Gelman and Shalizi, 2013; Gao et al., 2019). The gain in reducing bias in estimates and

increasing their precision was most evident with over-sampling (Tables 1, 2).

Our model could also provide weekly bycatch estimates which were largely unbiased except in the under-sampling scenario where a positive and negative bias remained at the beginning and end of a year respectively, even with 15 years of data (Figure 6). With under-sampling, few observed bycatch events can be collected by design because observers are very unlikely to be accepted on board by skippers. Weekly estimates were too high at the beginning of a year but too low at the end, but this somewhat canceled out at the year-level. There was still a slight overestimation bias resulting from our choice of a non-symmetric pattern for dolphin presence and a symmetric pattern for biased coverage: observing bycatch events at the end of a year was comparatively more difficult than at the beginning of a year because overlap between a non-nil coverage and dolphin presence was smaller at the end of year (Figure 1). These shortcomings illustrate that a model-based approach should
<table><caption>TABLE 2 | Statistical properties of estimates from the model-based approach.</caption><thead><tr><th>Method</th><th>Uncertainty</th><th>Data sampling</th><th>n<sub>years</sub></th><th>Bias (bycatch events)</th><th>Coverage (%)</th><th>Width of CI (%)</th><th>n<sub>obs</sub></th></tr></thead><tbody><tr><td>Model-based</td><td>Bayesian</td><td>Unbiased</td><td>3</td><td>3.0</td><td>100.0</td><td>91.1</td><td>45</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Unbiased</td><td>5</td><td>2.1</td><td>100.0</td><td>76.3</td><td>75</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Unbiased</td><td>10</td><td>1.1</td><td>100.0</td><td>59.1</td><td>150</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Unbiased</td><td>15</td><td>1.9</td><td>100.0</td><td>50.9</td><td>225</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Under-</td><td>3</td><td>10.0</td><td>100.0</td><td>164.6</td><td>15</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Under-</td><td>5</td><td>6.4</td><td>100.0</td><td>142.0</td><td>25</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Under-</td><td>10</td><td>8.3</td><td>100.0</td><td>112.9</td><td>50</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Under-</td><td>15</td><td>5.3</td><td>100.0</td><td>97.8</td><td>75</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Over-</td><td>3</td><td>7.4</td><td>100.0</td><td>53.2</td><td>63</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Over-</td><td>5</td><td>4.8</td><td>100.0</td><td>42.6</td><td>126</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Over-</td><td>10</td><td>3.5</td><td>100.0</td><td>32.6</td><td>630</td></tr><tr><td>Model-based</td><td>Bayesian</td><td>Over-</td><td>15</td><td>3.3</td><td>100.0</td><td>27.7</td><td>756</td></tr></tbody></table>

Several years of data were simulated a 100 times. Bias of point estimate (in number of bycatch events), coverage of (Bayesian 80%) credible interval (CI) and precision (as CI width relative to the point estimate) are reported. The last column indicates the average number of bycatch events ($n_{obs} = \frac{E}{\sum_{pt} Y_{pt}}$) that were recorded by onboard observers during data sampling.

be tailored to the context of the study, and we designed our simulations largely from our knowledge on the common dolphin in the Bay of Biscay. However, the framework of regularized multilevel regression with post-stratification is very flexible and we believe our proposed model has large potential for generality as it simply translates a decomposition of bycatch risk into a smooth time-varying and (unstructured) time-invariant effects. The model can easily be made more complex, data permitting, to accommodate spatial effects with, for example, a Besag-type prior (Sims et al., 2008; Morris et al., 2019).

Several important assumptions are structurally built into our model: in particular, a first order random was assumed for the mean function of the Gaussian Process prior, with no attempt to estimate from data the correlation parameter (e.g., using an AR(1) prior instead). The choice of a first order random walk was not aiming at uncovering the true data-mechanism: our aim were to reveal a temporal pattern in bycatch risk from sparse data using a flexible, yet parsimonious approach. This was particularly true in the under-sampling scenario where few bycatch events could be observed in any given year of simulated data. In the other scenarios, other choices than the first order random walk could be considered as more data are collected. We also assumed that the range parameter of the covariance function in the Gaussian Process prior for week effects was known and such that bycatch risk was temporally uncorrelated after 4 weeks. Fixing the range parameter is usually not recommended but was motivated by consideration of the data-to-parameter ratio, and computation convenience. Bycatch data are binary and can be sparse: these two features underscore how little information may be available. In this context, limiting the number of parameters to estimate can be justified on pragmatic consideration. The model we are proposing is parameter-rich, but some structure are assumed on these parameters in the form of the prior used. These priors represent choices from the analyst and

may be reconsidered and tested, data permitting. There was some evidence that bycatch risk was under-smoothed in the over-sampling scenario which resulted in an over-estimation of bycatch risk (Figure 5, rightmost panel). Model expansion is seamless with Stan (Gabry et al., 2019), and the above mentioned parameters could be estimated, rather than fixed, with adequate data. Despite somewhat arbitrary prior and modeling choices, our model provided more accurate estimates of bycatch numbers and bycatch risk in under- and over-sampling scenarios. This satisfactory predictive ability points to another important limitation.

Our model is phenomenological, i.e., it is agnostic of the causes behind the temporal variations in bycatch risk. Bycatch risk is the product of dolphin presence and bycatch probability given presence (the latter was constant in our simulations). The model only estimates this product of two probabilities and thus cannot disentangle them without other sources of data. This limitation seems inconsequential in our simulations for the aim of accurate estimation of the total number of bycatch events as interest lies in the effects of causes (how much bycatch?) rather than in the causes of effects (why bycatch occurred?). A straightforward model expansion (as pointed out by a reviewer) would be the consideration of *p* vessel-level covariates ($z_{1j},...,z_{pj}$) in Equation (1):

$$ \alpha_j \sim N \left( \sum_{b=1}^{p} (\xi_b \times z_{jb}), \sigma_{\text{vessel}} \right) \quad (7) $$

Candidate covariates such as vessel length or gear-attributes (e.g., mesh size) could be incorporated in the analysis to improve the exchangeability assumption on vessel-effects. An obvious covariate to consider for detecting self-selection of skippers into observer programme participation is to include whether a skipper
has ever accepted an observer, or the number of times it did so in the past: a negative regression coefficient could be interpreted as voluntary skippers having an intrinsically lower risk of bycatch. Including skipper-level covariates could reduce the between-skipper variance $\sigma_{vessel}^2$, and improve ultimately precision of bycatch estimates. Consideration of other distributions than the normal (e.g., a skew-normal, or a Student-t distribution with a fixed degree of freedom) would be straightforward with Stan but is probably worthwhile only with large enough amount of data for all practical purposes (McCulloch and Neuhaus, 2011).

An important assumption underlying accurate estimation is that the information on the total effort must also be accurate and available at the scale of weeks for post-stratification. This assumption is crucial to scale-up estimates from the (potentially biased) sample to the population, but it does not necessarily hold with fisheries effort as the latter is more often estimated rather than measured directly (Julian and Beeson, 1998; ICES, 2018, 2020b). Here we assumed that the total number of fishing operations (e.g., number of tows for trawls; Tremblay-Boyer and Berkenbusch, 2020) are available as auxiliary information for post-stratification. This assumption about the availability of disaggregated data stems from the explicit consideration of time as an important predictor of variations in bycatch risk. This assumption is necessary for using post-stratification to align the sample with the population targets but may be difficult to meet in practice. Currently, ICES WGBYC uses in its BRA a coarse, but admittedly comparable proxy across fisheries and countries to quantify fishing effort, namely days at sea (ICES, 2019). A day at sea is any continuous period of 24 h (or part thereof) during which a vessel is present within an area and absent from port (Anonymous, 2019a). Importantly, this definition is not at the level of a fishing operation, and effort thus quantified is already aggregated at a level above that at which bycatch data are collected. This coarsening of fisheries effort data is fundamentally a measurement problem, and one that modeling should not be expected to remedy easily. BRA uses an estimate of total fishing effort for the fisheries of concern in a specific region, together with some estimate of likely or possible bycatch rates that might apply for the species of concern, in order to evaluate whether or not the total bycatch in that area might be a conservation issue. A regularized multilevel regression model could be used to obtain estimates of bycatch rates to be used in BRA. Post-stratification could also be attempted using the coarse days at sea proxy for effort, and thus our framework could be adapted to match the requirements of ICES WGBYC.

Assuming that our framework were to be adopted to produce bycatch estimates, how would both fisheries and Non-Governmental Organizations (NGOs) react given the salience of bycatch as a policy issue in Europe? Such a prospective question inevitably entails some speculations (as with all “what-if” questions), but may nevertheless bring some insights as highlighted by a reviewer. Within Europe, the conservation reference currently available for assessing bycatch is that established under the Agreement on the Conservation of Small Cetaceans of the Baltic, North East Atlantic, Irish and North Seas. The agreement has the conservation objective to minimize

anthropogenic removals of harbor porpoises (and other small-sized cetaceans), and to restore and/or maintain population depletion to/at 80% or more of the carrying capacity in each assessment unit (ASCOBANS, 2000; ICES, 2020c). Methods for setting conservation reference points were agreed in March 2021 at the meeting of the Biodiversity Committee of the Olso-Paris Regional Sea Convention. This committee adopted the use of the Removals Limit Algorithm for harbor porpoises in the North Sea assessment unit and a modified Potential Biological Removal (Wade, 1998) for common dolphins in the North-East Atlantic (Genu et al.)<fnref n="3" />. Accurate bycatch estimates will be needed for assessment against these reference points. However, fisheries may challenge the accuracy of estimates precisely because they will result from a new statistical model. While a healthy skepticism is warranted, and model improvements are certainly possible, it must be kept in mind that our model only addresses the issue of having a correlation between observer coverage and bycatch risk, and does so with some assumptions. There would remain many biases to be addressed in bycatch data (Babcock and Pikitch, 2003), and many of them would be best addressed with a proper random allocation of professional observers to vessels (that is better design and better measurement). A purely model-based solution can be brittle (Sarewitz, 1999), and may lead to displacement of the problem of bycatch assessment to a never-ending problem of model improvement that would delay any corrective measures or decision (Rayner, 2012). Model-based estimates offer a pragmatic approach to the analysis of already collected data, but should not deflect from improving survey design where possible. Assuming that model-based estimates would be endorsed by a fishery industry, NGOs could challenge in court any reference point that is not zero for PETS, since by definition, it ought to be zero. The Habitats Directive requires strict protection and prohibits “all forms of deliberate capture or killing” (emphasis added) of all species listed on its Annex IV which includes all cetacean species. The Court of Justice of the European Union has consistently ruled that the adjective “deliberate” is to be understood in the sense of “conscious acceptance of consequences” (Trouwborst and Somsen, 2019): in other words, using knowingly a gear that may potentially catch a protected species contravenes the Habitat Directives. What will eventually play out remains to be seen, but strongly hinges on how polarized the bycatch issue is. As scientists, our duty remains to provide the best available evidence on bycatch and to outline all management actions and their consequences in light of this evidence (Pielke, 2007). Our model is unlikely to change bycatch management in France in the near term: both fisheries and NGOs are at loggerheads, vying for public and official support. They are building constituencies and advertising unyielding positions in diverse medias: we content that a legal confrontation at a national or supra-national level is extremely likely and probably being prepared. We nevertheless think our model, by making use of data already collected within the DCF framework and by encouraging further, ideally dedicated, monitoring; can be part of

<footnote n="3">³Genu, M., Gilles, A., Hammond, P., Macleod, K., Paillé, J., Paradinas, I. A., et al. (in preparation). Evaluating strategies for managing anthropogenic mortality on marine mammals: an R Implementation with the Package RLA.</footnote>
a messy solution to the wicked problem (Frame, 2008) of dolphin bycatch in the medium to long term, once the gavel hits and the dust settles.

## 5. CONCLUSION

We investigated with simulations the ability of multilevel regularized regression with post-stratification to estimate cetacean bycatch for observer programmes when coverage is correlated to bycatch risk. Our aims were to provide a first investigation on model-based estimates obtained from samples preferentially collected either during low- or high-bycatch risk periods. The unbiased sampling case is unrealistic (Babcock and Pikitch, 2003): biased sampling, either under-sampling or over-sampling (ICES, 2009), may be the general case. We considered both of these cases, under quite extreme scenarios whereby data collection was highly correlated with bycatch risk, resulting in either very few observed events with under-sampling, and a large number of observed events with over-sampling. In both cases, multilevel regularized regression with post-stratification was able to produce nearly unbiased bycatch estimates with as few as 5 observed events data. With only 1 year of data, precision was low, especially with under-sampling, and there was some estimation bias with over-sampling one. These results stemmed from the extreme scenarios we considered but illustrate nevertheless that a model cannot be expected to solve all the deficiencies of data collection and measurement. Good measurement is key for accurate estimation and our results actually re-emphasize the importance of design. However, they also show that a good data collection design and an adequate modeling framework are synergistic and allow to extract a lot of information for sparse data. Assuming a normal distribution for the bycatch estimates (which is not necessary as the posterior is available, but the following are back-of-the-envelope calculations to be used for deriving heuristics), a 80% Bayesian CI width divided by 2.5 gives an idea of the associated coefficient of variation: the model-based approach can yield a coefficient of variation of 50% with as few as 15 observed events if sampling is unbiased. With under-sampling, one would need 10 years of data (under our data simulation schemes) to obtain the same precision. This re-iterates the need to (i) have dedicated observer schemes, (ii) ensure adequate observer coverage and (iii) use a model-based approach tailored to extract as much information as possible from sparse data, as the first two points are very difficult to live up to in practice.

The key assumptions behind regularized multilevel regression with post-stratification in our simulations are that bycatch risk changes smoothly through time and that accurate data on the number of fishing operations at the same temporal scale are available (e.g., number of tows for trawls; Tremblay-Boyer and Berkenbusch, 2020). When both assumptions can be reasonably entertained, we showed how a model-based approach using recent methodological developments is attractive, irrespective of how data were collected. A further asset of the explicit consideration of a temporal scale is that it may help in pinpointing more precisely windows of heightened risk in order

to target adequate mitigation measures (e.g., spatio-temporal closures). The framework of multilevel modeling is very flexible and can accommodate spatial effects, etc., data permitting. Regularization will, in general, be needed to mitigate data sparsity and leverage partial pooling in order to obtain stable estimates of bycatch. Given the satisfactory performance of regularized multilevel regression with post-stratification in our simulations, we recommend further investigations using this technique to estimate bycatch rate and numbers from both representative or non-representative samples. The modeling choices we made (e.g., a first order random walk for the mean function, or fixing the range parameter in the covariance function of the Gaussian Process prior) are not prescriptive, and other choices of prior models for parameters should be investigated. Investigations should be tailored to the context, and modeling choices motivated by the latter: given the complexity of PETS bycatch, a one-size-fits-all solution is unlikely. A re-analysis of > 15 years of observer data on common dolphin bycatch in pair trawlers flying the French flag is currently underway (Rouby et al.)<fnref n="4" /> in order to obtain better bycatch estimates that could be further used to estimate conservation reference points in order to better manage this fishery in the long run (Cooke, 1999; Punt et al., 2021).

## DATA AVAILABILITY STATEMENT

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://gitlab.univ-lr.fr/mauthier/regularized_bycatch.

## AUTHOR CONTRIBUTIONS

MA: led the analyses, the conception, and writing of the paper. ER and KM: support in analyses, paper conception, and writing. All authors contributed to the article and approved the submitted version.

## FUNDING

ADERA and JNCC provided support in the form of salaries for MA and KM, respectively, but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## ACKNOWLEDGMENTS

We thank two reviewers from critical and constructive comments that help in improving the present manuscript. MA thanks Amaterasu for support when carrying the simulations.

<footnote n="4">⁴Rouby, E., Dubroca, L., Cloâtre, T., Demanêche, S., Genu, M., Macleod, K., et al. (in preparation). Estimating cetacean bycatch from non-representative samples (II): a Case Study on Common Dolphins in the Bay of Biscay.</footnote>
REFERENCES

Amandé, M. J., Chassot, E., Chavance, P., Murua, H., de Molina, A. D., and Bez, N. (2012). Precision in bycatch estimates: the case of Tuna Purse-Seine Fisheries in the Indian Ocean. *ICES J. Mar. Sci.* 69, 1501–1510. doi: 10.1093/icesjms/fss106

Anonymous (2016). Rapport Annuel sur la mise en œuvre du Règlement Européen (CE) 812/2004 Établissant les Mesures Relatives Aux Captures Accidentelles de Cétacés dans les Pêcheries. Rapport annuel réglementaire. Technical report, Direction des pêches maritimes et de l'aquaculture (DPMA).

Anonymous (2019a). Commission Implementing Decision (EU) 2019/910. Establishing the Multiannual Union Programme for the Collection and Management of Biological, Environmental, Technical and Socioeconomic Data in the Fisheries and Aquaculture Sectors. European Commission.

Anonymous (2019b). Regulation (EU) 2019/1241 of the European Parliament. The Council of 20 June 2019 on the Conservation of Fisheries Resources and the Protection of Marine Ecosystems Through Technical Measures.

ASCOBANS (2000). Annex O - report of the IWC-ASCOBANS Working Group on Harbour Porpoises. *J. Cetacean Res. Manage.* 2, 297–305.

Authier, M., Aubry, L., and Cam, E. (2017a). Wolf in sheep's clothing: model misspecification undermines tests of the neutral theory for life histories. *Ecol. Evol.* 7, 3348–3361. doi: 10.1002/ece3.2874

Authier, M., Saraux, C., and Péron, C. (2017b). Variable selection and accurate predictions in habitat modelling: a shrinkage approach. *Ecography* 40, 549–560. doi: 10.1111/ecog.01633

Avila, I. C., Kaschner, K., and Dormann, C. F. (2018). Current global risks to marine mammals: taking stock of the threats. *Biol. Conserv.* 221, 44–58. doi: 10.1016/j.biocon.2018.02.021

Babcock, E., and Pikitch, E. K. (2003). *How Much Observer Coverage is Enough to Adequately Estimate Bycatch?* Technical report, Pew Institute for Ocean Science.

Benoît, H., and Allard, J. (2009). Can the data from at-sea observer surveys be used to make general inference about catch composition and discards? *Can. J. Fish. Aquat. Sci.* 66, 2025–2039. doi: 10.1139/F09-116

Brownell, R. J., Reeves, R., Read, A., Smith, B., Thomas, P., Ralls, K., et al. (2019). Bycatch in gillnet fisheries threatens critically endangered small cetaceans and other aquatic megafauna. *Endangered Species Res.* 40, 285–296. doi: 10.3354/esr00994

Cam, E. (2012). "Each site has its own survival probability, but information is borrowed across sites to tell us about survival in each site": random effects models as means of borrowing strength in survival studies of wild vertebrates. *Anim. Conserv.* 15, 129–132. doi: 10.1111/j.1469-1795.2012.00533.x

Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., et al. (2017). Stan: a probabilistic programming language. *J. Stat. Softw.* 76, 1–32. doi: 10.18637/jss.v076.i01

Carretta, J. V., and Moore, J. E. (2014). *Recommendations for Pooling Annual Bycatch Estimates when Events are Rare*. Technical report. NOAA-TMNMS-SWFSC-528, NOAA. NOAA Technical Memorandum.

Carretta, J. V., Moore, J. E., and Forney, K. A. (2017). *Regression Tree and Ratio Estimates of Marine Mammal, Sea Turtle, and Seabird Bycatch in the California Drift Fishery: 1990-2015*. Technical Report NOAA-TM-NMFS-SWFSC-568, NOAA. NOAA Technical Memorandum.

Cloâtre, T. (2020). *Projet LICADO - Limitation des Captures Accidentelles de Dauphins Commun dans le Golfe de Gascogne. Bilan des Enquêtes auprès des Professionnels du Golfe de Gascogne* (Paris).

Cooke, J. G. (1999). Improvement of fishery-management advice through simulation testing of harvest algorithms. *ICES J. Mar. Sci.* 56, 797–810. doi: 10.1006/jmsc.1999.0552

Course, G. P., Pierre, J., and Howell, B. K. (2020). *What's in the Net? Using Camera Technology to Monitor, and Support Mitigation of, Wildlife Bycatch in Fisheries*. Technical report, Worldwide Wildlife Fund.

Dolman, S., Evans, P., Ritter, F., Simmonds, M., and Swabe, J. (2020). Implications of new technical measures regulation for cetacean bycatch in European waters. *Mar. Policy* 124:104320. doi: 10.1016/j.marpol.2020.104320

Downes, M., Gurrin, L. C., English, D. R., Pirkis, J., Currier, D., Spittal, M. J., et al. (2018). Multilevel regression and poststratification: a modeling approach to estimating population quantities from highly selected survey samples. *Am. J. Epidemiol.* 187, 1780–1790. doi: 10.1093/aje/kwy070

Fay, R., and Herriot, R. (1979). Estimates of income for small place: an application of james-stein procedures to census data. *J. Am. Stat. Assoc.* 74, 269–277. doi: 10.1080/01621459.1979.10482505

Forsberg, O. J. (2020). Polls and the US presidential election: real or fake? *Significance* 17, 6–7. doi: 10.1111/1740-9713.01437

Frame, B. (2008). 'Wicked,' 'Messy', and 'Clumsy': long-term frameworks for sustainability. *Environ. Plan C Govern. Policy* 26, 1113–1128. doi: 10.1068/c0790s

Gabry, J., Simpson, D., Vehtari, A., Betancourt, M., and Gelman, A. (2019). Visualization in Bayesian workflow. *J. R. Stat. Soc. Ser. A* 182, 389–402. doi: 10.1111/rssa.12378

Gao, Y., Kennedy, L., Simpson, D., and Gelman, A. (2019). *Improving Multilevel Regression and Poststratification with Structured Priors*. Technical report, University of Toronto & Columbia University.

Gelman, A. (2007). Struggles with survey weighting and regression modelling. *Stat. Sci.* 22, 153–164. doi: 10.1214/088342307000000203

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2014). *Bayesian Data Analysis, 3rd Edn.* Boca Raton, FL: CRC Press. doi: 10.1201/b16018

Gelman, A., Hill, J., and Vehtari, A. (2021). *Regression and Other Stories, 1st Edn.* Cambridge, MA: Cambridge University Press. doi: 10.1017/9781139161879

Gelman, A., and Shalizi, C. (2013). Philosophy and the practice of Bayesian statistics. *Brit. J. Math. Stat. Psychol.* 66, 8–38. doi: 10.1111/j.2044-8317.2011.02037.x

George, E. I., and Ročková, V. (2021). Comment: regularization via Bayesian penalty mixing. *Technometrics* 62, 438–442. doi: 10.1080/00401706.2020.1801258

Goldin, N., and Purse, B. V. (2016). Fast and flexible Bayesian species distribution modelling using gaussian processes. *Methods Ecol. Evol.* 7, 598–608. doi: 10.1111/2041-210X.12523

Gray, C. A., and Kennelly, S. J. (2018). Bycatches of endangered, threatened and protected species in marine fisheries. *Rev. Fish Biol. Fish.* 28, 521–541. doi: 10.1007/s11160-018-9520-7

Griffin, J., and Brown, P. (2017). Hierarchical shrinkage priors for regression models. *Bayesian Anal.* 12, 135–159. doi: 10.1214/15-BA990

Halsny, V. (2020). Nonresponse bias in inequality measurement: cross-country analysis using luxembourg income study surveys. *Soc. Sci.* Q. 101, 712–731. doi: 10.1111/ssqu.12762

Hansen, M. H., and Hurwitz, W. N. (1946). The problem of non-response in sample surveys. *J. Am. Stat. Assoc.* 41, 517–529. doi: 10.1080/01621459.1946.10501894

He, Y., Hodges, J., and Carlin, B. (2007). Re-considering the variance parametrization in multiple precision models. *Bayesian Anal.* 2, 529–556. doi: 10.1214/07-BA221

Hoffman, M. D., and Gelman, A. (2014). The no-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. *J. Mach. Learn. Res.* 15, 1593–1623.

ICES (2009). *Report of the Study Group for Bycatch of Protected Species (SGBYC)*. Technical report, Copenhagen. doi: 10.1093/ppar/19.1.22

ICES (2018). *Report from the Working Group on Bycatch of Protected Species, (WGBYC)*. Technical report, International Council for the Exploration of the Sea, Reykjavik.

ICES (2019). *Working Group on Bycatch of Protected Species (WGBYC)*. Technical report, International Council for the Exploration of the Sea, Faro.

ICES (2020a). *Bycatch of Protected and Potentially Vulnerable Marine Vertebrates - Review of National Reports under Council Regulation (EC) No. 812/2004 and Other Information*. Technical report.

ICES (2020b). *Report from the Working Group on Bycatch of Protected Species, (WGBYC)*. Technical report. International Council for the Exploration of the Sea.

ICES (2020c). *Workshop on Fisheries Emergency Measures to minimize BYCATCH of Short-Beaked Common Dolphins in the Bay of Biscay and Harbour Porpoise in the Baltic Sea (WKEMBYC)*. International Council for the Exploration of the Sea.

Jaramillo-Legorreta, A. M., Cardenas-Hinojosa, G., Nieto-García, E., Rojas-Bracho, L., Thomas, L., Ver Hoef, J. M., et al. (2019). Decline towards extinction
of Mexico's Vaquita Porpoise (*Phocoena sinus*). R. Soc. Open Sci. 6:190598.
doi: 10.1098/rsos.190598

Julian, F., and Beeson, M. (1998). Estimates of marine mammal, turtle, and
seabird mortality for two California Gillnet Fisheries: 1990-1995. *Fish. Bull.* 96,
271–284.

Kiewiet De Jonge, C. D., Langer, G., and Sinozich, S. (2018). Predicting state presidential election results using national tracking polls and multilevel regression with poststratification (MRP). *Public Opin. Q.* 82, 419–446.
doi: 10.1093/poq/nfy023

Komoroske, L. M., and Lewison, R. L. (2015). Addressing fisheries bycatch in a
changing world. *Front. Mar. Sci.* 2:83. doi: 10.3389/fmars.2015.00083

Laran, S., Authier, M., Blanck, A., Dorémus, G., Falchetto, H., Monestiez, P., et al. (2017). Seasonal distribution and abundance of cetaceans within French waters-part II: the Bay of Biscay and the English channel. *Deep Sea Res. II* 141, 31–40. doi: 10.1016/j.dsr2.2016.12.012

Lax, J., and Phillips, J. (2009). How should we estimate public opinion in the states? *Am. J. Polit. Sci.* 53, 107–121. doi: 10.1111/j.1540-5907.2008.00360.x

Lennert, C., Kruse, S., Beeson, M., and Barlow, J. (1994). *Estimates of incidental marine Mammal Bycatch in California Gillnet Fisheries for July Through December, 1990*. Report of the International Whaling Commission.

Martin, R., Pearmain, E. J., Burfield, I. J., Small, C., Phillips, R. A., Yates, O., et al. (2019). Threats to seabirds: a global assessment. *Biol. Conserv.* 237, 525–537. doi: 10.1016/j.biocon.2019.06.033

Martin, S. L., Stohs, S. M., and Moore, J. E. (2015). Bayesian inference and assessment for rare-event bycatch in marine fisheries: a drift gillnet fishery case study. *Ecol. Appl.* 25, 416–429. doi: 10.1890/14-0059.1

McCulloch, C. E., and Neuhaus, J. M. (2011). Misspecifying the shape of a random effects distribution: why getting it wrong may not matter. *Stat. Sci.* 26, 388–402.
doi: 10.1214/11-STS361

Monnahan, C. C., Thorson, J. T., and Branch, T. A. (2017). Faster estimation of bayesian models in ecology using Hamiltonian Monte Carlo. *Methods Ecol. Evol.* 8, 339–348. doi: 10.1111/2041-210X.12681

Morris, M., Wheeler-Martin, K., Simpson, D., Mooney, S. J., Gelman, A., and DiMaggio, C. (2019). Bayesian hierarchical spatial models: implementing the Besag York Mollié model in Stan. *Spat Spatio Temp Epidemiol.* 31:100301. doi: 10.1016/j.sste.2019.100301

National Marine Fisheries Service (2004). *Evaluating Bycatch: a National Approach to Standardized Bycatch Monitoring Programs*. Technical Report. U.S. Department of Commerce, NOAA.

National Marine Fisheries Service (2005). *Revisions to Guidelines for Assessing Marine Mammal Stocks*. Technical report, NOAA.

National Marine Fisheries Service (2016). *Guidelines for Preparing Stock Assessment Reports Pursuant to the 1994 Amendments to the MMPA*. Technical report, NOAA.

Neal, R. (1998). "Regression and classification using Gaussian process priors," in *Bayesian Statistics*, eds J. M. Bernardo, J. O. Berger, A. P. Dawid, and A. F. M. Smith (Oxford: Oxford University Press), 475-501.

Northridge, S., Kingston, A., and Thomas, L. (2019). *Annual Report on the Implementation of Council Regulation (EC) No 812/2004 During 2018*. Technical report, University of St Andrews.

Pacoureau, N., Rigby, C. L., Kyne, P. M., Sherley, R. B., Winker, H., Carlson, J. K., et al. (2021). Half a century of global decline in oceanic sharks and rays. *Nature* 589, 567–571. doi: 10.1038/s41586-020-03173-9

Palka, D. L., and Rossman, M. C. (2001). *Bycatch Estimates of Coastal Bottlenose Dolphin (Tursiops truncatus)* in US Mid-Atlantic Gillnet Fisheries for 1996 to 2000. Technical Report. Northeast Fisheries Science Center.

Peltier, H., Authier, M., Caurant, F., Dabin, W., Daniel, P., Dars, C., et al. (2021). In the wrong place at the wrong time: identifying spatiotemporal co-occurrence of bycaught common dolphins and fisheries in the Bay of Biscay (NE Atlantic) from 2010 to 2019. *Front. Mar. Sci.* 8:617342. doi: 10.3389/fmars.2021.617342

Peltier, H., Authier, M., Deaville, R., Dabin, W., Jepson, P., van Canneyt, O., et al. (2016). Small cetacean bycatch as estimated from stranding schemes: the common dolphin case in the Northeast Atlantic. *Environ. Sci. Policy* 63, 7–18. doi: 10.1016/j.envsci.2016.05.004

Pérez, M., Perrichi, L. R., and Ramírez, I. (2017). The scaled beta2 distribution as a robust prior for scales. *Bayesian Anal.* 12, 615–637. doi: 10.1214/16-BA1015

Pielke, R. A. Jr. (2007). *The Honest Broker - Making Sense of Science in Policy and Politics, 1st Edn.* Cambridge, MA: Cambridge University Press.

Plummer, M., Best, N., Cowles, K., and Vines, K. (2006). CODA: convergence diagnosis and output analysis for MCMC. *R News* 6, 7–11.

Punt, A. E., Siple, M., Francis, T. B., Hammond, P. S., Heinemann, D., Long, K. J., et al. (2021). Can we manage marine mammal bycatch effectively in low-data environments? *J. Appl. Ecol.* 58, 596–607. doi: 10.1111/1365-2664.13816

R Core Team (2020). *R: A Language and Environment for Statistical Computing*. Vienna: R Foundation for Statistical Computing.

Rayner, S. (2012). Uncomfortable knowledge: the social construction of ignorance in science and environmental policy discourses. *Econ. Soc.* 41, 107–125. doi: 10.1080/03085147.2011.637335

Rhem, J., Kilian, C., Rovira, P., Shield, K. D., and Manthey, J. (2020). The elusiveness of representativeness in general population surveys for alcohol. *Drug Alcohol* 40, 161–165. doi: 10.1111/dar.13148

Sala, A., Konrad, C., and Doerner, H. (eds.). (2019). *Review of the Implementation of the EU Regulation on the Incidental Catches of Cetaceans (STECF-19-07)*. Publications Office of the European Union, Luxembourg.

Sarewitz, D. (1999). How science makes environmental controversies worse. *Environ. Sci. Policy* 7, 385–403. doi: 10.1016/j.envsci.2004.06.001

Sims, M., Cox, T., and Lewison, R. (2008). Modeling spatial patterns in fisheries bycatch: improving bycatch maps to aid fisheries management. *Ecol. Appl.* 18, 649–661. doi: 10.1890/07-0685.1

Stan Development Team (2020). *RStan: the R interface to Stan*. R package version 2.21.22. New York, NY.

St. John, F., Keane, A., Jones, J., and Milner-Guiland, E. (2014). Robust study design is as important on the social as it is on the ecological side of applied ecological research. *J. Appl. Ecol.* 51, 1479–1485. doi: 10.1111/1365-2664.12352

Stransky, C., and Sala, A. (2019). *Scientific, Technical and Economic Committee for Fisheries (STECF) - Revision of the EU-MAP and Work Plan Template (STECF-19-12)*. Technical report, Publications Office of the European Union, Luxembourg.

Tregenza, N. J. C., Berrow, S. D., Hammond, P. S., and Leaper, R. (1997). Harbour porpoise (*Phocoena phocoena L.*) by-catch in set gillnets in the Celtic Sea. *ICES J. Mar. Sci.* 54, 896–904. doi: 10.1006/jmsc.1996.0212

Tremblay-Boyer, L., and Berkenbusch, K. (2020). *Characterisation of Marine Mammal Interactions with Fisheries & Bycatch Mitigation*. Report for project INT2019-03 prepared for Department of Conservation. Technical report, Dragonfly Data Science.

Trouwborst, A., and Somsen, H. (2019). Domestic cats (*Felis catus*) and European nature conservation law-applying the EU birds and habitats directives to a significant but neglected threat to wildlife. *J. Environ. Law* 32, 391–415. doi: 10.1093/jel/eqz035

Ulrich, C., and Doerner, H. (eds.). (2021). *Scientific, Technical and Economic Committee for Fisheries (STECF) - 66th Plenary Report (PLEN-21-01)*. Publications Office of the European Union, Luxembourg.

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P. (2019). *Rank-Normalization, Folding, and Localization: an Improved *R* for Assessing Convergence of McMC*. Technical report, Helsinki Institute for Information Technology, Department of Computer Science, Aalto University.

Wade, P. R. (1998). Calculating limits to the total allowable human-caused mortality of cetaceans and pinnipeds. *Mar. Mammal Sci.* 14, 1–37. doi: 10.1111/j.1748-7692.1998.tb00688.x

Wang, W., Rothschild, D., Goel, S., and Gelman, A. (2015). Forecasting elections with non-representative polls. *Int. J. Forecast.* 31, 980–991. doi: 10.1016/j.ijforecast.2014.06.001

Zahorski, A. (2020). *Multilevel Regression With Post-Stratification for the National Level VIBER/Street Poll on the 2020 Presidential Election in Belarus*. Technical report. Uladzimir Karatkevich National University of Belarus.

Zhang, X., Holt, J. B., Lu, H., Wheaton, A. G., Ford, E. S., Greenland, K. J., et al. (2014). Multilevel regression and poststratification for small-area estimation of population health outcomes: a case study of chronic obstructive pulmonary disease prevalence using the behavioural risk factor surveillance system. *Am. J. Epidemiol.* 179, 1025–1033. doi: 10.1093/aje/kwu018
**Conflict of Interest:** MA is employed by the commercial company ADERA which did not play any role in this study beyond that of employer.

The remaining authors declare that the research was conducted in the absence of
any commercial or financial relationships that could be construed as a potential
conflict of interest.

**Publisher's Note:** All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in

this article, or claim that may be made by its manufacturer, is not guaranteed or
endorsed by the publisher.

Copyright © 2021 Authier, Rouby and Macleod. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
<figure><img src="image_227.png" /></figure>

# Estimating Bycatch From Non-representative Samples (II): A Case Study of Pair Trawlers and Common Dolphins in the Bay of Biscay

Etienne Rouby<fnref n="1" />,²,³, Laurent Dubroca⁴, Thomas Cloâtre⁵, Sebastien Demanèche⁵, Mathieu Genu<fnref n="1" />,³, Kelly Macleod⁶, Hélène Peltier³,⁷, Vincent Ridoux<fnref n="1" />,²,³ and Matthieu Authier<fnref n="1" />,³,⁷*

OPEN ACCESS

**Edited by:**
Jeff Moore,
Southwest Fisheries Science Center
(NOAA), United States

**Reviewed by:**
Stephen Stohs,
Southwest Fisheries Science Center
(NOAA), United States
Tomo Eguchi,
National Marine Fisheries Service
(NOAA), United States
Doug Sigourney,
Integrated Statistics, United States

*Correspondence:
Matthieu Authier
matthieu.authier@univ-lr.fr

†Present address:
Kelly Macleod,
HIDef Aerial Surveying Limited,
Lillyhall, United Kingdom

**Specialty section:**
This article was submitted to
Marine Megafauna,
a section of the journal
Frontiers in Marine Science

Received: 15 October 2021
Accepted: 02 December 2021
Published: 13 January 2022

**Citation:**
Rouby E, Dubroca L, Cloâtre T, Demanèche S, Genu M, Macleod K, Peltier H, Ridoux V and Authier M (2022) Estimating Bycatch From Non-representative Samples (II): A Case Study of Pair Trawlers and Common Dolphins in the Bay of Biscay. Front. Mar. Sci. 8:795942. doi: 10.3389/fmars.2021.795942

Marine megafauna plays an important functional role in marine ecosystems as top predators but are threatened by a wide range of anthropogenic activities. Bycatch, the incidental capture of non-targeted species in commercial and recreational fisheries, is of particular concern for small cetacean species, such as dolphins and porpoises. In the North-East Atlantic, common dolphin (*Delphinus delphis*, Linné 1758) bycatch has been increasing and associated with large numbers of animals stranding during winter on the French Atlantic seashore since at least 2017. However, uncertainties around the true magnitude of common dolphin bycatch and the fisheries involved have led to delays in the implementation of mitigation measures. Current data collection on dolphin bycatch in France is with non-dedicated observers deployed on vessels for the purpose of national fisheries sampling programmes. These data cannot be assumed representative of the whole fisheries' bycatch events. This feature makes it difficult to use classic ratio estimators since they require a truly randomised sample of the fishery by dedicated observers. We applied a newly developed approach, regularised multilevel regression with post-stratification, to estimate total bycatch from unrepresentative samples and total fishing effort. The latter is needed for post-stratification and the former is analysed in a Bayesian framework with multilevel regression to regularise and better predict bycatch risk. We estimated the number of bycaught dolphins for each week and 10 International Council for the Exploration of the Sea (ICES) divisions from 2004 to 2020 by estimating jointly bycatch risk, haul duration, and the number of hauls per days at sea (DaS). Bycatch risk in pair trawlers flying the French flag was the highest in winter 2017 and 2019 and was associated with the longest haul durations. ICES divisions 8.a and 8.b (shelf part of the Bay of Biscay) were estimated to have the highest common dolphin bycatch. Our results were consistent with independent estimates of common dolphin bycatch from strandings. Our method shows cases how non-representative observer data can nevertheless be analysed to estimate fishing duration, bycatch risk and, ultimately, the

<footnote n="1">¹ Université de La Rochelle, La Rochelle, France, ² UMR7372 Centre d'Etudes Biologiques de Chizé (CEBC), Villiers En Bois, France, ³ Observatoire PELAGIS, UMS3462, CNRS-La Rochelle Université, La Rochelle, France, ⁴ Ifremer, Laboratoire Ressources Halieutiques de Port en Bessin, Port-en-Bessin-Huppain, France, ⁵ Ifremer, Laboratoire de Biologie Halieutique de Brest, Brest, France, ⁶ Joint Nature Conservation Committee, Aberdeen, United Kingdom, ⁷ ADERA, Pessac, France</footnote>
number of bycaught dolphins. These weekly-estimates improve upon current knowledge of the nature of common dolphin bycatch and can be used to inform management and policy decisions at a finer spatio-temporal scale than has been possible to date. Our results suggest that limiting haul duration, especially in winter, could serve as an effective mitigation strategy.

**Keywords:** additional mortality, anthropogenic activities, modelling, non-representative samples, conservation, small cetaceans, fisheries, post-stratification

# 1. INTRODUCTION

Over the last 50 years, the conservation status of cetaceans has been deteriorating (Brownell et al., 2019). Over 80 species of cetaceans occur worldwide and bycatch, the non-intentional capture or killing of non-target species in commercial or recreational fisheries (Hall, 1996; Davies et al., 2009), remains a threat, especially to small-sized species (Scarfiff, 1977; Read et al., 2006; Avila et al., 2018; Anderson et al., 2020). Success stories in small cetacean conservation are the exception rather than the rule (e.g., Bessesen, 2018). Both Rogan et al. (2021) and Bearzi and Reeves (2021) opined of institutional failures to conserve cetaceans in European Waters in spite of current legislation (for example, the Habitats Directive, the Marine Strategy Framework Directive) or regional agreements such as the Agreement on the Conservation of Small Cetaceans of the Baltic and North Seas (ASCOBANS, see **Table 1** for acronyms; ICES, 2020c). Over 20 species of small cetaceans have been registered in the North-East Atlantic, with roughly half of which occurring regularly (Course, 2021). Because of their slow life histories and their limited potential rates of increase, small cetaceans are particularly at risk of decline when anthropogenic activities induce additional mortality on populations (Read, 2008). Anthropogenic activities and their cumulative impacts can take a heavy toll on populations. Common species may disappear, such as short-beaked common dolphins (*Delphinus delphis*, hereafter called common dolphins) in the Adriatic Sea (Bearzi and Reeves, 2021), or are under many threats, e.g., in the Bay of Biscay (García-Baron et al., 2019; Murphy et al., 2021).

In 2013, the common dolphin's conservation status in the European Marine Atlantic, as assessed under Article 17 of the Habitats Directive, was “Unfavourable–Inadequate” because of fishery bycatch (Murphy et al., 2021). Common dolphin bycatch in the Bay of Biscay, in particular, has attracted a lot of media coverage since 2017 in international outlets<fnref n="1" /> and motivated (with bycatch of Harbour porpoise *Phocoena phocoena* in the Baltic Sea) a special request of Non-Governmental Organisations to the European Commission in 2019. The International Council for the Exploration of the Sea (ICES) advised in 2020, for the common dolphin in the Bay of Biscay, a combination of temporal closures of all métiers (i.e the combination of gear, target species, and fishing area) of concern and application of pingers on pair trawlers to mitigate bycatch outside of the period of closure (ICES, 2020b). Temporal closures, restricted to winter months in

<table><caption>**TABLE 1** | List of acronyms.</caption><thead><tr><td>Acronym</td><td>Meaning</td></tr></thead><tbody><tr><td>ASCOBANS</td><td>Agreement on the Conservation of Small Cetaceans of the Baltic, North East Atlantic, Irish and North Seas</td></tr><tr><td>DCF</td><td>Data collection framework</td></tr><tr><td>DPMA</td><td>“Direction des pêches maritimes et de l’aquaculture”</td></tr><tr><td>GNS</td><td>Gillnetters</td></tr><tr><td>GTR</td><td>Gill trammel netters</td></tr><tr><td>ICES</td><td>International Council for the Exploration of the Sea</td></tr><tr><td>Ifremer</td><td>Institut Français de Recherche pour l’Exploitation de la Mer</td></tr><tr><td>ObsMer</td><td>Observation des captures en Mer (French national observer scheme for monitoring fisheries)</td></tr><tr><td>PBR</td><td>Potential Biological Removal</td></tr><tr><td>PTM</td><td>Pair trawlers</td></tr><tr><td>PTB</td><td>Bottom pair trawlers</td></tr><tr><td>VAST</td><td>Vector-Autoregressive Spatio-Temporal</td></tr><tr><td>WGBYC</td><td>ICES Working Group on Bycatch of Protected Species</td></tr><tr><td>WKEMBYC</td><td>ICES Workshop on fisheries Emergency Measures to minimise BYCatch of short-beaked common dolphins in the Bay of Biscay and harbor porpoise in the Baltic Sea</td></tr></tbody></table>

which strandings of common dolphins with evidence of bycatch have increased in recent years (ICES, 2020d), could have been implemented as emergency measures under the provisions of the Common Fisheries Policy. For 2021, France instead required the mandatory use of acoustic repulsive devices (pingers) on all pair trawlers flying the French Flag (code métier Pair trawlers and hereafter referred to as PTM) operating in the Bay of Biscay<fnref n="2" />, a technical mitigation measure whose efficiency was found wanting (Ulrich and Doerner, 2021). This decision against the advice of ICES was motivated by a lack of knowledge on common dolphins, including its abundance at the level of the whole North-East Atlantic (the currently recognised management unit: Murphy et al., 2013) and the extent of bycatch. The issue of managing uncomfortable knowledge through interpretation of scientific uncertainty can be raised (Schweder, 2000; Rayner, 2012); yet it should not eclipse that there are genuine difficulties in estimating accurately the true magnitude and the extent of bycatch of small cetaceans (Moore et al., 2021).

Several types of fishing gear are known to cause cetacean bycatch: drift nets, set gill, trammel nets, both pair and single midwater trawls, and some demersal trawls (Rogan and Mackey,

<footnote n="1">¹https://www.nytimes.com/2019/05/02/world/europe/france-dolphins-fishing.html</footnote>

<footnote n="2">²https://www.mer.gouv.fr/protecting-cetaceans-annick-giradin-presents-7-commitments-made-french-state-fishermen-and</footnote>
2007; Fernández-Contreras et al., 2010; Peltier et al., 2016). Accurate quantification of bycatch rates by fishing gears or métiers remains a challenging endeavour (Babcock et al., 2003; ICES, 2019). Traditionally, bycatch data are collected by onboard observers monitoring fishing operations and recording the unwanted catch of non-commercial species (Course, 2021). Ratio estimators, based on the number of observed hauls with bycatch over the total number of monitored hauls, are used (Alverson et al., 1994; page 18) but are plagued by large uncertainties due to low coverage and the usual small number of hauls with small cetacean bycatch (Babcock et al., 2003; Authier et al., 2021; Course, 2021). It may also happen that some bycatch events may not be reported by non-dedicated observers since they may drive observations for other purposes than report bycatch (e.g., commercial discards or stock assessments). A critical assumption behind the use of such ratio-estimators is that of a representative sample: this assumption is difficult to sustain unless monitoring is dedicated to marine mammals, and allocation of observers to fishing vessels is truly randomised (that is, not at the discretion of skippers). Even if we are willing to assume representative sampling, if coverage is low, the main challenge remains to extrapolate from sample to the whole fisheries. In France, monitoring of cetacean bycatch in fisheries is non-dedicated (Cornou et al., 2018), and the collected data are described as non-representative of the bycatch events, preventing the use of ratio-estimators (Anonymous, 2016; page 24).

This non-dedicated nature and the sparseness of the bycatch data complicates the use of state-of-the-art spatio-temporal models such as Vector-Autoregressive Spatio-Temporal (VAST) (Thorson, 2019). This framework accommodates density-dependence, spatial and temporal scales to estimate biomass or abundance or presence of a species (Thorson et al., 2015). Spatio-temporal models are also used to model the co-occurrence of commercial and bycaught species, allowing the estimate of bycatch risk with time-varying spatial effects (Ward et al., 2015). These types of model-based approaches methodologies allow modelling spatial and temporal auto-correlation through the use of Gaussian process priors. It is difficult to transfer a priori the same model-based structure to analyse small cetacean bycatch. Models such as VAST capitalise on the availability of catch data that are collected as part of fisheries monitoring. In contrast, bycatch monitoring is not as developed or efficiently enforced in many fisheries in Europe (ICES, 2019, 2020a; Sala et al., 2019), and bycatch data are typical of low quality and unrepresentative (Authier et al., 2021). In Europe, fisheries monitoring is carried out under the Data Collection Framework (DCF) but "remains not well-suited for the dedicated monitoring of rare and protected bycatch in high-risk fisheries since its main focus is the statistically-sound random sampling of all commercial fisheries" (Ulrich and Doerner, 2021). Because of these data quality issues, Authier et al. (2021) conducted a simulation study to gauge the potential of investigating recent methods for the analysis of non-representative samples (for a recent example of a model-based approach to estimate bycatch, refer to Luck et al., 2020) in the context of small cetacean bycatch: they concluded the potential of regularised multilevel regression with post-stratification to infer

more accurately bycatch rates (although uncertainties remained large). The approach of Authier et al. (2021) also makes use of Gaussian process priors but does not necessarily assume that a large dataset has been collected.

We analysed historical bycatch monitoring data collected by onboard observers (from 2004 to 2020) on PTM, a métier historically associated with high levels of dolphin bycatch in the Bay of Biscay (ICES, 2019; Murphy et al., 2021). Leveraging recent modelling developments (see companion article; Authier et al., 2021), we jointly estimated bycatch risk, haul duration, and number of hauls per days at sea (DaS) from an updated and revised observer dataset on common dolphin bycatch. The modelling procedure accounts for the sparseness of the bycatch incident dataset and the low observer coverage through constraints. This type of constraint (which can be viewed as some sort of penalisation) is also called regularisation. We used structured priors, such as Gaussian processes, to achieve regularisation and leverage the within-year information at the weekly scale (inducing correlation between some weeks). Structured priors allow inducing some spatial- or temporal-dependency between so called random-effects whereas unstructured priors do not induce such dependency (but both assume exchangeability). Importantly, we used this model-based approach to disaggregate bycatch risk at the level of calendar weeks in order to document within-year variations. Estimates were summed over a whole year to investigate between-year variations in the number of bycaught dolphins. We compared these model-based estimates with strandings, both within- and between-years. Finally, we concluded with recommendations on conservation and mitigation.

# 2. MATERIALS AND METHODS

## 2.1. Materials

### 2.1.1. Study Division

The study area (Figure 1) encompasses 10 ICES divisions within area 27: it includes the Bay of Biscay, the English Channel, and part of the Celtic seas. These zones are associated with submesoscale and mesoscale oceanographic processes, such as eddies and upwelling, that enhance ecosystem productivity and result in high availability of fishes, including commercial species (e.g., European seabass *Dicentrarchus labrax*, Sardine *Sardina pilchardus* or Anchovy *Engraulis encrasicolus*). Each division can roughly be classified as oceanic or neritic: divisions 7.d, 7.e, 7.f, 7.g, 7.h, 8.a, 8.b, and 8.c are related to neritic ecosystems while divisions 7.j, and 8.d are related to oceanic ecosystems.

### 2.1.2. Data Sources

Two main sources of data were used. The first dataset, called ObsMer<fnref n="3" /> ("Observation des captures en Mer"), is collected as part of an onboard observer program set up within the Data Collection Framework of the Common Fisheries Policy. The ObsMer program is carried out by Ifremer ("Institut Français de Recherche pour l'Exploitation de la Mer"), under the supervision of the Directorate of Fisheries and Aquaculture ("Direction

<footnote n="3">³https://sih.ifremer.fr/Ressources/ObsMer</footnote>
<figure><img src="image_228.png" /><figcaption>**FIGURE 1** | Study area in the North-East Atlantic ocean, with ICES divisions overlayed.</figcaption></figure>

des pêches maritimes et de l'aquaculture," DPMA). ObsMer observers' primary duty is to register the length and weight composition of catches. Still, they have to report any bycatch event if they witness such events. ObsMer data on PTM cover 4,484 hauls between 2004 and 2021, of which 82 were associated with a bycatch event of at least 1 and up to 50 common dolphins. ObsMer provides, among other information, the geographic position, timing, and duration of hauls. Although ObsMer is aiming at a coverage of 10 and 5% of fishing effort for (level-3 métier) PTM for vessels of more than 15 m and less than 15 m, respectively, these figures are rarely, if ever, reached in practice: accepting onboard observers remains entirely at the discretion of skippers. The effort is quite low overall, ranging from 0 to 11% of Days at Sea (DaS) (Table 2). A DaS is any continuous period of 24 h (or part thereof) during which a vessel is present within an area and absent from the port (Anonymous, 2019). The number of observed hauls with at least one bycatch record is very small because the yearly percentage of observed hauls with a bycatch event never exceeded 4.5% and was 0 in nearly half of the surveyed years. ObsMer data on pair-trawlers are an unrepresentative sample of hauls, largely because allowing an observer remains largely at the discretion of skippers (Babcock et al., 2003; Benoît and Allard, 2009).

The second dataset provides monthly estimates of total fishing effort in each division. This dataset is generated from the algorithm SACROIS developed by Ifremer and integrates data from Vessel Monitoring System, log-books, and landing statistics (for boats longer than 18 m from January 1, 2004, and longer than 15 m from January 1, 2005; Système d'Information Halieutique, 2017). SACROIS aims at (1) correcting errors that could exists in the integrated dataset due to recording or collecting errors and

(2) reconstitute métiers during the fishing trip as they are not recorded in logbooks or fish market data (Cornou et al., 2018). The SACROIS dataset provides the best available estimates of total effort, in DaS, between 2004 and 2020 (Table 2). There are also refusals from skippers due to administrative and security reasons. Skippers must file an application for authorisation to embark observers and even if they decide to file, the authorisation may be declined due to security reasons (e.g., not enough room or rails not high enough).

These two datasets are complementary for our purposes: ObsMer provides micro-level data on marine mammal bycatch at the resolution of hauls. From these data, bycatch risk may be estimated (Luck et al., 2020). Fishing trips effort data, on the other hand, are macro-level: they provide spatialised effort data at the scale of a whole fishing fleet. These population-level data on effort allows the post-stratification of bycatch risk estimated from observer data to obtain the number of bycaught dolphins (Authier et al., 2021). Descriptive statistics of both datasets are displayed in Table 2. Used in tandem, both datasets allow using regularised multilevel regression with post-stratification to estimate cetacean bycatch from non-representative samples (Authier et al., 2021).

Estimates were finally compared to strandings along the French Atlantic seaboard. The French Stranding Network, founded in the 1970s, is dedicated to the monitoring of marine mammal strandings along the shores of France (mainland and overseas). Around 400 trained volunteers are currently taking an active part in the network. These volunteers make the complete coverage of French coastlines possible. Standardised training of volunteers by permanent Observatoire Pelagis staff, which takes place two times a year, ensure the homogeneity, comparability, and standardisation of data collection procedures in the field. Observatoire Pelagis is mandated by the French Ministry of Ecology to train and deliver authorisation to handle carcasses of marine mammals (which are all protected species under national law). It also collates the data and analyse it to inform on the status of marine mammal populations. Stranding data for the period 2004-2020 were used. Only common dolphins found with lesions diagnostic of bycatch in fishing gear were considered (Kuiken, 1994) as well as those stranded during multiple stranding events, or "unusual mortality events" related to lesions diagnostic of bycatch. Multiple stranding events were defined as high numbers of strandings occurring in a restricted area with a common cause of death. The threshold was defined at 30 cetaceans over 10 consecutive days recorded along a maximal distance of 200 km in the Bay of Biscay, and 10 individuals per 10 days per 200 km of coastline along the coast of the western Channel (Peltier et al., 2014). Reverse drift modelling uses a deterministic drift model developed by Meteo France (Peltier et al., 2012) to reconstruct the trajectory of every stranded common dolphin from its stranding location to its likely area of death at sea. The number of dead stranded animals in each cell is then corrected by the cell-specific probability of being stranded (Peltier et al., 2016). These probabilities were estimated by numerical experiment in which the drift of carcasses in the study area was simulated in order to assess with which frequency they would reach a coastline (Peltier and Ridoux, 2015).