Introduction
The success of conservation efforts relies foremost on our ability to quantify variability in the abundance and distribution of wildlife at a variety of spatial and temporal scales. Effective species conservation and management also requires information on species distribution patterns, but this can be particularly challenging for highly mobile and cryptic species that may be subject to multiple anthropogenic stressors across international boundaries. Understanding species–habitat relationships can improve the assessment of population size, trends, and distribution by explicitly allowing high‐resolution data on habitats to inform abundance estimation and the identification of protected areas. Worldwide, the marine environment is under increasing pressure from human activities (Halpern et al. , Maxwell et al. ). Top predators, such as marine mammals, can serve as prime sentinels of multitrophic marine ecosystem changes since they integrate ecological variations across large spatial and long temporal scales (Moore , Kiszka et al. ).
The European Union (EU) has recognized increasing human pressures (e.g., intended and unintended harvest, habitat degradation, underwater noise, and climate change) to the marine environment and adopted policies, such as the Marine Strategy Framework Directive (MSFD; EU‐COM ), to manage human impacts and to conserve biodiversity (Bertram and Rehdanz ). However, information on the spatio‐temporal distribution of marine mammals throughout the European seas is often lacking. The harbor porpoise (Phocoena phocoena) is among the most common cetaceans in European Atlantic shelf waters (Hammond et al. , ) and as a top predator considered an important indicator species. It is distributed in coastal and shelf waters and is particularly exposed to human impacts that can lead to disturbance, injury, or death (e.g., Siebert et al. , Vinther and Larsen , Das et al. , Dähne et al. , Thompson et al. , Dyndo et al. ). In EU waters, the species has a high conservation status under several conventions; for example, it is listed in Annexes II and IV of the EU Habitats Directive (EU‐COM ) and on the OSPAR list of threatened and/or declining species and habitats of the Northeast Atlantic (OSPAR ). Harbor porpoises within European waters show distinct variation in their distribution on seasonal, inter‐annual, and decadal time scales (Gilles et al. , Scheidat et al. , Hammond et al. , Peschko et al. ); however, drivers of these variations have not yet been fully identified, in part because ecological linkages between porpoises, their environment and their prey species are dynamic and poorly understood. Species–habitat models (or species distribution modeling, SDM) can help shed light on ecological relationships and/or provide a basis for hypothesis‐driven studies on trophic linkages.
The ecological theory of SDM assumes that species distributions are largely influenced by environmental variables (Austin ). Although migration, predator avoidance, and social interactions influence the behavior of cetaceans, distribution patterns are most strongly influenced by foraging on patchy prey (Palacios et al. , , Redfern et al. ). Most marine mammals aggregate in spatially constrained areas, often termed “hotspots,” where prey availability and foraging efficiency are high. Prey availability in the ocean, in turn, is largely determined by physical forcing mechanisms, such as fronts and eddies that aggregate primary and secondary producers. For harbor porpoises, both static variables (e.g., water depth) and dynamic variables (e.g., currents, tide‐related variables, chlorophyll) have been shown to explain spatial distribution patterns (Marubini et al. , Embling et al. , Gilles et al. , Booth et al. , IJsseldijk et al. ). The incorporation of habitat variables into cetacean distribution models allows the estimation of abundance and animal densities at finer spatial and temporal scales than conventional, design‐based line‐transect analyses (see reviews in Redfern et al. , Gregr et al. ); therefore, such model‐based densities can improve fine‐scale management needs.
Previously, the only robust estimates of harbor porpoise density for the entire North Sea were derived from the two SCANS (Small Cetacean Abundance in the North Sea) surveys, conducted in 1994 and 2005 (Hammond et al. , ). While there has been considerable regional survey effort on national scales (e.g., Gilles et al. , Haelters et al. , Scheidat et al. , Geelhoed et al. , Peschko et al. ), these data sets have been analyzed separately and provide a spatially and temporally patchy picture of porpoise distribution in the south‐eastern North Sea. The aim of this study was to combine the available comparable data sets to develop seasonal habitat‐based density models for harbor porpoises in the central and southern North Sea. The data sets were collected over 9 yr (2005–2013) and during three seasons by means of dedicated aerial surveys for harbor porpoises that used standardized line‐transect survey methods and incorporated correction factors for missed animals on the transect line (Hiby and Lovell , Buckland et al. ). We developed generalized additive models of porpoise density, incorporating near real‐time remote sensing data to identify dynamic mesoscale oceanographic features that are often observed to be hotspots for marine mammals (Scales et al. ). The results of our study are especially relevant for marine spatial planning, which requires accurate fine‐scale maps of species distribution to assess risks of increasing human activities at sea.
Materials and Methods
Study area
Our study area covered 411,000 km2 in the central North Sea and Southern Bight ranging from 3.4° W to 9.8° E and from 50.9° N to 58.4° N (Fig. ). The North Sea is a continental shelf sea, characterized by shallow water depths (30–200 m), except for the Norwegian Channel and the Skagerrak (up to 600 m) on its eastern margin (Fig. ). The highest productivity occurs in coastal regions and in areas such as the Dogger Bank, which is the largest sandbank in the central North Sea (18,000 km2), where high levels of phytoplankton production occur year‐round (Brockmann et al. ). The North Sea is rich with frontal patterns, and up to 10 oceanic fronts were distinguished by Belkin et al. (). Tidal mixing, river plume, and upwelling fronts dominate these areas (Krause et al. ). A high number of sea surface temperature (SST) fronts were detected (Belkin et al. ). Other front‐like structures are the boundaries between the North Sea water and the Baltic outflow from the Skagerrak, as well as between the inflow of Atlantic water from the North and from the English Channel (Otto et al. ).
Study area in the North Sea. (A) Water depth and distribution of sandeel grounds; (B) aerial survey strata in Denmark (DK), Germany (DE), the Netherlands (NL), and Belgium (BE); (C) survey strata covered during the SCANS II (ship and aerial) and the aerial Dogger Bank surveys. The dashed line indicates the Exclusive Economic Zones (EEZ).
Survey data and field methodology
We based our density models on a large number of dedicated aerial visual surveys conducted within the respective national waters of Denmark, Germany, the Netherlands, Belgium and in adjacent areas (UK) between 2005–2013 (Figs. and ). All study areas, except for the smaller Belgian EEZ (see Haelters et al. ), were geographically stratified (Fig. ) and a parallel or zigzag transect layout was used to provide equal coverage probability in each survey block. All these systematic surveys were part of national monitoring programs conducted throughout most of the year, although at varying frequencies that led to uneven regional coverage. To increase the spatial coverage of the central North Sea, we included several blocks of shipboard and aerial survey from the SCANS II survey conducted during summer 2005 (Hammond et al. , Fig. ).
Overview of available overall monthly (top panel) and yearly (bottom panel) survey data in the period 2005–2013.
Distance‐sampling line‐transect methods (Buckland et al. ) were consistently followed on all the visual surveys and the same experienced observers participated in most of the surveys. Detailed descriptions of aerial survey field methods are provided elsewhere (Scheidat et al. , Gilles et al. , Hammond et al. ). All aerial surveys used the same platform, a twin engine, high‐wing aircraft outfitted with two bubble windows for unobstructed downward and lateral viewing. Surveys were flown at an altitude of 600 ft (183 m) and groundspeeds of 90–100 kts (165–185 km/h). Two observers scanned the water surface through left and right bubble windows. A third person, the data recorder, entered all reported data directly into a laptop with real‐time GPS input using a customized program (except in Belgium where the observers recorded sighting information manually). The aircraft's position was stored every 2 s. The start and end positions of the transect lines and the exact sighting positions were recorded. Surveys were only conducted in Beaufort sea states 0–3 and visibilities > 5 km. Sea state, turbidity, cloud coverage and, for each observer side, glare details (angle obscured and intensity of glare) and an overall subjective assessment of detection conditions (good, moderate, or poor) were recorded at the beginning of each transect and whenever conditions changed. Good conditions normally required sea states ≤ 2, clear water and no strong glare; conditions changed to moderate when sea state deteriorated or strong glares affected downward view. Total effective strip widths (esw; both sides of the plane), taking account of detection probability less than 1 on the transect line (commonly known as g(0)), were estimated for good and moderate conditions using the racetrack data collection method (Hiby and Lovell , Hiby ). This provided a correction for missed animals within varying sighting conditions that could also differ between observer sides (see Scheidat et al. for details). Here, total esw was estimated to be 170 m (CV = 0.23) under good conditions and 67 m (CV = 0.24) under moderate conditions, incorporating g(0) values of 0.405 and 0.164, respectively.
Shipboard cetacean survey data gathered during the SCANS II survey in blocks U and V were used to increase sampling coverage in the central North Sea. The shipboard surveys used double platform line‐transect survey methods, described in Hammond et al. (). In order to correct for animals missed on the transect line, we applied a g(0) of 0.216 (CV = 0.16) for the SCANS II ship survey blocks (taken from Hammond et al. ), but re‐estimated the stratum (ship) specific esw for blocks U and V.
Data preparation
Survey data collected in poor conditions (as defined for aerial surveys) and in Beaufort sea states higher than 2 for shipboard surveys were not included. Survey effort was very low and patchy during winter; therefore, data collected during December–February were excluded from further analysis. The three remaining seasons were defined according to the meteorological start in north temperate zones (i.e., Mar–May = spring; Jun.–Aug. = summer; Sep.–Nov. = fall); this division also corresponds well to temperature cycles in the North Sea (Nauw et al. ). In addition, this division also captures sensitive periods in the annual reproduction cycle of harbor porpoises in the North Sea, where the birth and mating period occurs in summer (Lockyer ).
The SCANS II ship data were already prepared and divided into segments of continuous effort over approximately 30 min, resulting in average segment length of 6.8 km (SD = 2.8; median = 7.34 km). The aerial survey transects were divided into continuous‐effort segments with a target segment length of 10 km, to capture the characteristics of the environmental data and to reduce the number of segments with zero sightings. Leftover segments of less than 2 km were ignored (i.e., 359 segments with 24 sightings), as commonly done (e.g., Marubini et al. , Booth et al. , Mannocci et al. ). Since this affected < 2% of the data, the omission of a short segment at the end of a transect would not be expected to introduce a bias.
Predictor variables
The candidate set of predictor variables (Table ) was selected based on ecological knowledge and previous studies on the distribution of harbor porpoises. Due to large data gaps, remotely sensed chlorophyll‐a concentration was excluded from the list of initially considered variables. The set included the projected spatial covariates x and y, bathymetry and topography related predictors, closest distance to coast, distance to sandeel (Ammodytes spp.) fishing grounds and, as dynamic predictors, sea surface temperature (SST), the spatial and temporal standard deviation (SD) of SST (as proxies for fronts), and day length to capture seasonal variations (Table ).
List of candidate predictors and their abbreviation as used throughout the text. All predictors were included in the model without any transformationPredictor | Abbreviation | Unit | Description |
x ETRS coordinate | x | m | ETRS Lambert Azimuthal, EPSG 3035 |
y ETRS coordinate | y | m | |
Water depth | Depth | m | Bathymetrie |
Slope | Slope | Degree | Slope of seabed; slope and aspect were computed using “terrain” in R package raster 2.3–24 (Hijmans 2015a) |
Aspect | Aspect | Degree | Azimuthal direction of slope |
Distance to sandeel ground | Dist to sandeel | km | Jensen et al. , closest distances were retrieved using rgeos 0.3–8 (Bivand and Rundel ) in R |
Distance to coast | Dist to coast | km | Closest distance to European mainland or islands; rgeos 0.3–8 |
Sea surface temperature | SST | °C | Daily image at 1 km pixel resolution; blended high‐resolution MUR SST (JPL ) |
SST‐SD‐Time | SST‐SD‐Time | °C | SD of SST for 8‐d period ending on each survey day |
SST‐SD‐Space5 | SST‐SD‐Space5 | °C | SST variability in 11 × 11 km box |
SST‐SD‐Space10 | SST‐SD‐Space10 | °C | SST variability in 21 × 21 km box |
SST‐SD‐Space20 | SST‐SD‐Space20 | °C | SST variability in 41 × 41 km box |
Day length | Dayl | h | Photoperiod for given latitude and date; R package geosphere 1.3–13 (Hijmans 2015b) |
The sandeel is an important prey species not only for harbor porpoises (Santos and Pierce ) but also for piscivorous fish that are also preyed upon by porpoises. A recent study in the southern North Sea showed that sandeel was among the “big four” of harbor porpoise diet, besides gobies, gadoids, and clupeids (Leopold and Meesters ). In our model, data on the spatial distribution of sandeel foraging grounds were based on information from fisheries and navigation data provided by fishermen (Jensen et al. ). A distance of zero was assigned to all effort segments lying within these sandeel grounds.
Remotely sensed SST was derived on a daily basis at 0.011 degrees (about 1 km) equal‐angle horizontal resolution from a level 4 (L4; gap‐free gridded) SST image, downloaded from NOAA's ERDDAP data server (Simons ). Since oceanic features such as fronts often result in steep SST gradients, we included the SD of SST as a proxy for frontal regions. A spatial SD of SST (SST‐SD‐Space) was calculated for a square box with a radius of 5, 10, and 20 pixels in x and y direction centered on the pixel containing the segment midpoint. To capture temporal variability at a given point (as a proxy for frontal movement), we also included an 8‐d temporal measure of variation (SST‐SD‐Time).
Collinearity between predictors was addressed by estimating the variance inflation factor (VIF) using the R package HH 3.1‐15 (Heiberger ). VIF for a single explanatory variable was obtained using the r‐squared value of the regression of that variable against all other explanatory variables. Large VIF values indicate collinearity (Dormann et al. ). All candidate predictors yielded VIF scores below 4, suggesting that collinearity was not a major concern in our data, with highest VIFs for the spatial covariates x (2.92) and y (3.81) as well as distance to coast (2.95), slope (2.75), and water depth (2.67).
Model framework
Generalized additive models (GAMs; Hastie and Tibshirani ) have previously been used to model nonlinear relationships between cetaceans and habitat covariates (e.g., Becker et al. , , Gilles et al. , Forney et al. , ). GAMs were fitted in R 3.1.2 (R Core Team ) using the package mgcv 1.8.‐3 (Wood , ). Restricted maximum likelihood (REML) with automatic term selection (Marra and Wood ) was used for smoothing parameter estimation. We used default thin‐plate regression splines in all models except for the tensor product smooths of spatial interaction terms, where the default cubic regression splines were used.
The model was fit to the full data set (i.e., 2005–2013, all months from March‐Nov.). The number of porpoise groups encountered per segment was defined as the response variable. Models assuming a quasi‐Poisson, negative binomial, and Tweedie distribution for the response variable were initially considered. Based on inspection of diagnostic plots of model residuals and quantiles, the negative binomial distribution (‘nb()’ function in mgcv) was ultimately selected for the final models. The natural logarithm of the effective area searched (in km2) was included as an offset to account for both varying segment lengths and varying detection probabilities based on recorded sighting conditions during the surveys. Sighting condition specific estimates of esw were applied for each segment with good vs. moderate subjective sightings conditions, as defined above for aerial surveys. The effective area searched was then calculated as the segment length multiplied by esw (incl. g(0)). Since these estimates are specific to sighting conditions and platform, any differences in sighting rates/detection probabilities are fully accounted for.
As part of model diagnostics we also tested for spatial autocorrelation (SAC) in the model residuals using a correlogram to assess data independence (Dormann ).
Model selection
Model selection was performed in a two‐step process that first identified candidate predictors to include in the model using goodness of fit criteria (e.g., AIC, REML, and visual inspection of modeled spatial patterns compared to sighting data) and then selected a final model from the candidate models based on predictive performance (cross‐validation across years).
Step 1: Identify candidate models based on goodness of fit measures
- Candidate models were initially identified using AIC/REML criteria within mgcv. Models with and without a smooth of spatial covariates x and y were considered separately and with each of the SST‐SD‐Space predictors (5, 10, and 20) since these were highly correlated.
- Nonsignificant predictors (α = 0.05) and predictors where the confidence band of the partial residuals included zero throughout the model's predictor range (i.e., the predictor might not be influential) were subsequently excluded and models were re‐fit prior to the cross‐validation process.
- The prediction plots for the remaining candidate models were inspected and the average squared prediction error (ASPE; Hastie and Tibshirani ) was calculated across all segments to assess predictive performance on the full data set. Models were excluded from consideration if mean porpoise density predictions throughout the North Sea study area were excessively large and biologically infeasible, as this was considered an indication of poor model convergence or mis‐parameterization.
Step 2: Select the best candidate model based on cross‐validation
In a second step, we performed cross‐validation to identify the candidate model with the greatest predictive power on novel data sets, using a pseudo‐jackknife process similar to that described in Becker et al. (). Each candidate model was fit to subsets of the available survey data that either excluded 2009, 2010, 2011, 2012, or 2013, respectively. These years are characterized by high overall effort and good spatial coverage of the study area. Earlier years were not used as novel data sets for cross‐validation because they either had limited spatial coverage (2006–2008; Appendix S1: Fig. S1), or they represented the only coverage of the northern part of our study area and were considered essential to the model (2005). During the re‐fitting, all predictors that were included in the candidate models for the full data set were included. Model predictions for the excluded (“novel”) year were evaluated using ASPE. The model that minimized the sum of ASPE across all the above individual‐year cross‐validation steps was selected as the final “best” model.
Prediction
Predictions were made on a spatial grid of habitat and static covariates at a resolution of 5 × 5 km. Given the finer spatial resolution of the prediction grid, we also made predictions on a 10 × 10 km grid to evaluate potential bias. We found that the predictions for all three seasons were almost identical, but the coarser 10 × 10 km scale is less useful to managers. We limited the grid to an area of about 411,000 km2 that included all covered transects (Fig. ), to avoid predicting outside of the range of covariates used in model fitting. Porpoise group densities were predicted on a daily basis for each survey period in each year resulting in 1864 daily predictions. To allow comparison of seasonal trends in porpoise density, we subsequently averaged daily group densities within each season (spring, summer, and fall) and multiplied by the mean observed seasonal group size (Table ).
Summary of 2005–2013 survey data used for model development: effective survey effort (km in good or moderate conditions for the aerial surveys, and with sea state ≤2 for the ship surveys) as well as number of sightings of harbor porpoise (Hp) groups and mean group size are shown. Seasonal mean group sizes are based on the combined sighting data for each season (see text for definitions)Survey | Effort (km) | Hp sightings | Hp individuals | Mean group size (SD) | Source |
Belgium (2008–2013) | 9207 | 1192 | 1442 | 1.21 | 1, 2 |
Dogger Bank (2011 and 2013) | 12,246 | 1125 | 1708 | 1.52 | 3, 4 |
Denmark, south. North Sea (2011–2013) | 2561 | 227 | 291 | 1.28 | 4 |
Germany (2005–2013) | 87,263 | 8720 | 10,604 | 1.22 | 5, 6, 7, 8, 9 |
The Netherlands (2008–2013) | 38,013 | 2718 | 3244 | 1.19 | 10, 11 |
SCANS II aerial surveys 2005 (blocks B, H, L and Y) | 4248 | 225 | 277 | 1.23 | 12 |
SCANS II ship surveys 2005 (blocks U and V) | 3119 | 149 | 232 | 1.56 | 12 |
Sum | 156,630 | 14,356 | 17,798 | 1.24 (0.47) | |
Spring | 59,486 | 6461 | 7466 | 1.16 (0.37) | |
Summer | 71,430 | 6797 | 8842 | 1.30 (0.51) | |
Fall | 25,714 | 1098 | 1490 | 1.36 (0.60) |
Notes:
SD, standard deviation. Sources are: 1, Haelters et al. (); 2, Haelters et al. (); 3, Gilles et al. (2012a); 4, Geelhoed et al. (); 4, Hansen (); 5, Gilles et al. (); 6, Gilles et al. (); 7, Gilles et al. (2014a); 8, Gilles et al. (2014b); 9, Dähne et al. (); 10, Scheidat et al. (); 11, Geelhoed et al. (); 12, Hammond et al. ().
Model‐based abundance and estimation of uncertainty
Model‐based abundance was calculated as the product of the predicted animal density in each cell and the “at‐sea” area of the cell (in km2), summed over all grid cells to obtain an overall seasonal model‐based abundance estimate for the entire study area. Sources of uncertainty in our model‐based abundance estimates include variations in mean group size (CV1), combined model specification and encounter rate (CV2), and detection parameters (CV3 & CV4). Each of these components was estimated separately. Uncertainty in the estimated mean group size (CV1) was calculated from the sightings in each season (Table ). Following the approach in Forney et al. (), a jackknife procedure was used to estimate model uncertainty and encounter rate variation (CV2). For the jackknife, the survey days were randomly divided into 10 sets, and one set (comprising 10% of the survey days) was withheld for each of 10 jackknife iterations. Model fitting and prediction grid calculations were performed for each jackknife iteration as described above for the full data set, implying that uncertainty associated with the GAM coefficients is incorporated. Model uncertainty was then estimated on a cell‐by‐cell basis as the SD of the jackknife model predictions. Uncertainty in the detection parameters was estimated based on the standard error of esw and g(0), respectively, for aerial (CV3) and shipboard surveys (CV4; only in summer). To retrieve an overall coefficient of variation (CVall) and to calculate log‐normal confidence intervals (CI) for the seasonal model‐based abundance estimates, the different components of uncertainty were then combined following Goodman ():[Image Omitted. See PDF]
Eq. assumes independence between the different sources of uncertainty, which is justifiable in our case since there are no common parameters or shared data subsets for the different sources of uncertainty.
Results
Joined survey data
A total of 251 survey (effort) days conducted during 60 months within the study period 2005–2013 were included in this study. We aggregated 156,630 km of on‐effort survey data with 14,356 sightings of harbor porpoise groups to build the models (for details on survey blocks, transect lines, and sighting locations see references in Table and Fig. ).
A total of 16,673 effort segments were included in the full model, with a median segment length of 9.88 km. Of these, 6305 segments were assigned to spring, 7657 to summer and 2711 to fall (Fig. ). Data were available for the south‐eastern North Sea in spring and fall, whereas a larger portion of the North Sea was covered during summer due to the inclusion of the two Dogger Bank surveys and the large‐scale SCANS II survey that took place in summer (Fig. ). However, the northern (i.e., north of 56° N) and southwestern regions of the study area were rather poorly covered.
Seasonal coverage of transect segments in 2005–2013, for spring (Mar.–May), summer (Jun.–Aug.), and fall (Sep.–Nov.). Effort segments are shown in gray, sighting positions in red.
Model selection
The first step of model selection yielded a total of six candidate models (Table ). These models either included an isotropic bivariate function of spatial location (x and y) or a three‐dimensional tensor product and either distance to coast or depth or both. During the second step of model selection, model 4 was selected as the “best” model because it minimized the sum of ASPE across all individual‐year cross‐validation steps (Table ) and had the best model diagnostics and goodness of fit measures (i.e., lowest REML/AIC). This best model explained 24.7% of the deviance using a negative binomial error distribution and a three‐dimensional tensor product smooth of location and SST, a smooth of day length, distance to sandeel grounds, distance to coast, average water depth, SST‐SD‐Space20 and SST‐SD‐Time (ordered according to decreasing Chi‐square scores; Fig. ). This model also significantly outperformed a simple (null) model assuming a uniform distribution of porpoises. Slope dropped out during model selection, since the confidence band included zero throughout the covariate's range. The ratio of the observed to model‐predicted density estimates for the entire study area was 0.97. No significant autocorrelation in model residuals could be detected within the first 10 lags.
Candidate models fitted to the full suite of survey data, using the no. porpoise sightings as response variable, and associated goodness of fit measures. Theta = value of adjustment parameter identified from the negative binomial distribution. Offset = log(effective area searched in km2)Model | Theta | Model formula | REML | AIC (ΔAIC) | % Dev. | ASPE | Sum ASPE |
1 | 0.80 | s(x,y) + s(SST) + te(SST, dayl) + s(SST‐SD‐Time) + s(SST‐SD‐Space20) + s(depth) + s(dist to coast) + s(dist to sandeel) + offset | 19164.6 | 38175.7 (130.4) | 23.7 | 2.62 | 14.06 |
2 | 0.79 | s(x,y) + s(SST) + te(SST, dayl) + s(SST‐SD‐Time) + s(SST‐SD‐Space20) + s(dist to coast) + s(dist to sandeel) + offset | 19185.6 | 38225.9 (180.6) | 23.3 | 2.62 | 14.06 |
3 | 0.79 | s(x,y) + s(SST) + te(SST, dayl) + s(SST‐SD‐Time) + s(SST‐SD‐Space20) + s(depth) + s(dist to sandeel) + offset | 19169.5 | 38191.6 (146.3) | 23.5 | 2.62 | 14.00 |
4 | 0.83 | te(x,y,SST) + s(dayl) + s(SST‐SD‐Time) + s(SST‐SD‐Space20) + s(depth) + s(dist to coast) + s(dist to sandeel) + offset | 19107.7 | 38045.3 (0) | 24.7 | 2.60 | 13.84 |
5 | 0.82 | te(x,y,SST) + s(dayl) + s(SST‐SD‐Time) + s(SST‐SD‐Space20) + s(dist to coast) + s(dist to sandeel) + offset | 19124.1 | 38083.1 (37.8) | 24.4 | 2.60 | 13.86 |
6 | 0.81 | te(x,y,SST) + s(dayl) + s(SST‐SD‐Time) + s(SST‐SD‐Space20) + s(depth) + s(dist to sandeel) + offset | 19128.0 | 38117.6 (72.3) | 24.2 | 2.61 | 13.85 |
Notes:
REML, Restricted maximum likelihood; AIC, Akaike's information criterion; %Dev, percentage of deviance explained by the model. ASPE (average squared prediction error) was calculated across all segments within the entire study area. Sum ASPE = sum of ASPE across all individual‐year cross‐validation steps.
Functional plots of environmental variables relative to harbor porpoise density as indicated by the estimated smooth functions for the selected covariates in the best model. Plots of 1‐dim smooths are shown, whereas the 3‐dim tensor product smooth of location and SST (te(x, y, SST)) cannot be displayed. Estimated degrees of freedom (edf) for nonlinear fits are provided in parentheses on the y‐axes. Hatch marks on the x‐axes show sample values and range of samples. The shaded areas (2× standard error bands) denote the 95% Bayesian confidence intervals (CI). Note for interpretation that some CIs tend to be very large at the higher edges of the observed covariate values, where sampling was limited. In general, the range of partial residual values on the y‐axis provides a visual cue of how important a variable might be in concert with all the selected predictors.
Day length and the tensor product smooth of location and SST proved to be important predictors and good proxies for “season,” allowing dynamic predictions in both space and time. The inclusion of a three‐dimensional tensor product that includes spatial components (such as coordinates x and y) in conjunction with an ecologically sound variable is a common approach to capture local fluctuations of the respective ecological variable (e.g., Augustin et al. , Williams et al. ). Densities generally increased with day length, with highest densities predicted when day length exceeded 14.5 h during the months of June through August. Highest porpoise density occurred at 150 km offshore and at depths between 25 and 40 m (Fig. ). Furthermore, harbor porpoise densities increased with SST‐SD‐Time, that is, higher probability for SST‐fronts, and decreased with distance to sandeel grounds (Fig. ).
Seasonal predictions
Our spatial predictions of porpoise density across the central and southeastern North Sea are consistent with previously described seasonal patterns of porpoises in the respective national waters (e.g., Gilles et al. , , Haelters et al. , , Scheidat et al. , Geelhoed et al. ). The density map for the spring season showed major hotspots of harbor porpoises in the southern and southeastern part of the North Sea, mainly inshore close to the Belgian and Dutch coasts extending toward the German coast off the East Frisian Islands (Fig. ). The model predicted high densities in the area of the Sylt Outer Reef in the German North Sea as well as north off the coast of Jutland in Denmark. Another hotspot in spring was identified at the Dogger Bank and northwest of this large sandbank. No survey data were available in spring for this region, and hence the model predictions should be seen as spatial and temporal extrapolations. Standard deviations for both the Dogger Bank and northern Jutland areas were generally higher than those of other regions where survey data were available. The higher model uncertainties in these regions also resulted in larger 90% CIs (Fig. ). Density predictions for these areas should thus be viewed with caution.
Predicted harbor porpoise densities in the North Sea in spring (Mar.–May). Upper panel: The overlaid contours are associated jackknife standard deviations (SD), whereas the black and white dashed boundary depicts the sampling coverage in spring (concave hull of survey segments; also see Fig. ). Lower panel: Lower and upper lognormal 90% confidence intervals (Lower 90% CI and Upper 90% CI) for the seasonal density based on the jackknife samples.
In contrast to spring, harbor porpoise hotspots in summer shifted toward offshore and western areas, and a large hotspot was present off the German and Danish west coast that extended toward the Dogger Bank (Fig. ). The hotspot off the East and West Frisian Islands identified in spring was still present in summer, although smaller in size.
Predicted harbor porpoise densities in the North Sea in summer (Jun.–Aug.) Upper panel: The overlaid contours are associated jackknife standard deviations (SD). Lower panel: Lower and upper lognormal 90% confidence intervals (Lower 90% CI and Upper 90% CI) for the seasonal density based on the jackknife samples.
The predicted density in fall was lower than in the other seasons. The distribution was spatially heterogeneous and areas with higher densities were predicted north‐west of the Dogger Bank and off the German and Danish west coasts (Fig. ).
Predicted harbor porpoise densities in the North Sea in fall (Sep.–Nov.). Upper panel: The overlaid contours are associated jackknife standard deviations (SD), whereas the black and white dashed boundary depicts the sampling coverage in fall (concave hull of survey segments; also see Fig. ). Lower panel: Lower and upper lognormal 90% confidence intervals (Lower 90% CI and Upper 90% CI) for the seasonal density based on the jackknife samples.
Model‐based abundance
Seasonal model‐based predictions of abundance were similar in spring and summer; but the predicted abundance in fall was about a third lower than the summer/spring abundances (Table ). The difference among seasons was, however, not significant.
Model‐based abundance (N) estimate for harbor porpoises in the North Sea study area (as depicted in Figs. 5–7; area size = 411,000 km2). CVall (N) includes weighted measures of all sources of variance (model specification, mean group size, and detection probabilities)Parameter | Spring | Summer | Fall |
Abundance (N) | 372,167 | 361,146 | 228,913 |
CV (model) | 0.006 | 0.001 | 0.003 |
CV (group size) | 0.01 | 0.01 | 0.02 |
CV (weighted mean esw) | 0.18 | 0.12 | 0.19 |
CV (ship g(0)) | – | 0.16 | – |
CVall (N) | 0.18 | 0.20 | 0.19 |
Lower 95% CI | 260,658 | 243,827 | 159,264 |
Upper 95% CI | 531,380 | 534,913 | 329,022 |
Density | 0.91 | 0.88 | 0.56 |
Note:
CI, Confidence interval.
The same relationship was found when estimating model‐based abundance for the equivalent area covered in all three seasons (i.e., the area coverage in spring, see Fig. ; size = 119,240 km2), resulting in estimated abundances of 129,922 (95% CI: 91,556–184,365) in spring, 118,322 (95% CI: 80,259–174,437) in summer and 72,105 individuals in fall (95% CI: 49,850–104,295).
The sources of uncertainty in the abundance estimation differ substantially. The highest uncertainty in the abundance estimate is caused by high CVs for esw and g(0), whereas the uncertainty caused by group size estimation and model specification are relatively low (Table ).
We also compared the model‐based abundance estimates for all three seasons between two grid resolutions (5 × 5 and 10 × 10 km) and found that they were identical, consistent with what is expected from spatial models (Miller et al. ).
For further quality assurance we compared previously published design‐based standard line‐transect (LT) abundance estimates within individual national and SCANS II strata to model‐based estimates derived from the daily density predictions for the same time period in each stratum (Fig. ). The confidence limits for the model‐based abundance estimates, as predicted by the best model for certain spatio‐temporal subsets, overlap with the 95% CIs of LT estimate for all examined survey strata. The apparent large difference for SCANS II block U is attributable to the greater uncertainty in both estimates for this stratum (Fig. ).
Model‐based abundance estimates (model) for harbor porpoises in comparison to those derived from design‐based line‐transect (LT) surveys, plotted in chronological order (first by source/country then by month/year). Error bars represent 95% confidence intervals. See Fig. for abbreviations and locations of strata. Sources of LT estimates: (1) Hammond et al. , (2) Haelters et al. , (3) Haelters et al. , (4) Scheidat et al. , (5) Geelhoed et al. , (6) Gilles et al. 2012a, (7) Geelhoed et al. , (8) Dähne et al. , (9) Gilles et al. 2014a, (10) Gilles et al. 2012b, (11) Gilles et al. 2014b.
Discussion
Our habitat‐based density model provides fine‐scale information on the seasonal density and distribution of harbor porpoises in the North Sea. We built the model using an unprecedented data set that included systematic survey data collected in a large part of the North Sea over a 9‐yr period. By analyzing these data jointly we were able to significantly improve our understanding of which habitats and environmental conditions influence the distribution and abundance of harbor porpoises in the North Sea and to identify geographic regions of high porpoise density in different seasons. The novel, daily predictions of habitat‐based density allow the maximum flexibility to meet a variety of scales for management, as the daily predictions of porpoise density can be averaged over different time scales (e.g., weeks, months, seasons, years) depending on management needs; along with estimates of temporal variation and model uncertainty.
Since our surveys were repeated throughout the study period we were able to evaluate the predictive power of our habitat models (Guisan and Zimmermann ) on “novel” data sets collected in the same region but during another year, which revealed that our models were robust against missing years of data.
The applied methods are transferable to any other species and ecosystem, and by building on the experience and workflows from cetacean studies conducted in the north Pacific, we set an important milestone on how to standardize methods, with regard to model selection and validation, across habitats.
Previously, a broad‐scale prediction of harbor porpoise distribution in the complete North Sea was based on the SCANS II survey in July 2005 (Hammond et al. ). SCANS II indicated an increased porpoise density south of 56° N, approximately twice as much as estimated for the previous survey in 1994 (SCANS, Hammond et al. ). As the estimated total abundance did not change between 1994 and 2005, it was assumed that this was due to a southward shift in porpoise distribution. By combining the large‐scale SCANS II survey with the more frequent, fine‐scale national and Dogger Bank surveys, we were able to develop seasonal maps that improve our understanding of spatial and temporal variability in porpoise density and distribution in the North Sea. By shifting our focus from the national perspective to a more international perspective it becomes clear that similar seasonal inshore/offshore distribution patterns exist in adjacent areas.
What is striking is the comparable low density in fall. On one hand, this could be a sampling artifact as effort has been low and varied in spatially small areas (Figs. and ). As the survey data collected in fall were mostly based on national monitoring schemes, offshore areas were not covered and the predicted densities are at unconfirmed low levels. On the other hand, it is also possible that the movement pattern of the North Sea porpoises indeed leads to a—yet undocumented—northern (or southern) shift out of our study area. Unfortunately, only few publications exist on the distribution of harbor porpoises in offshore waters. The only satellite telemetry study with harbor porpoises in the area is Sveegaard et al. (), in which 64 porpoises were tagged in the Skagerrak and Inner Danish waters between 1997 and 2007. This study showed pronounced seasonal variations in the usage of several high‐density areas. However, most tagged individuals did not move into the central and southern North Sea, and none moved south of 54° N and into the German and Dutch EEZ (Sveegaard et al. ). In November 2011, an acoustic survey was conducted in the central North Sea across the Dogger Bank. Porpoises were detected acoustically throughout the study area; however, detection rates were higher over the UK region of the Bank than either the Dutch or German regions (MCR ). These findings provide support for the predicted hotspot west of the Dogger Bank in fall. However, we have very little means to assess the fall distribution of porpoises thoroughly and any future management action and research plan should address this essential data gap.
The role of environmental variables on porpoise density
This study is one of the first that included a measure of prey occurrence, namely sandeel, as a covariate in the porpoise density model. However, rather than using actual sandeel abundance, the distance to sandeel habitats was used. This approach acknowledges that sandeels are mobile and may move outside the known grounds. The sandeel is an important prey species for harbor porpoises and other top predators in the North Sea (e.g., Santos and Pierce , Leopold and Meesters ). Sandeels are forage fish, described as major food‐web energy conveyers (Christensen et al. ) and as such are important high‐energy mid‐trophic prey for piscivorous fish, such as gadoids, that are also porpoise prey.
An underlying hypothesis of our modeling approach and the selected candidate predictors is that harbor porpoises respond to spatial gradients of oceanographic variables. These relationships may vary across study areas, as underlying ecological processes and ranges of available predictor values differ. In this study, we focused on remotely sensed SST and temporal and spatial variation therein as the main dynamic predictors, because spatial gradients of SST indicate the presence of frontal systems where prey may congregate (Franks ). In addition, SST fronts are typically co‐located with fronts in other water properties such as salinity, density, and chlorophyll (Belkin et al. ). Indeed, studies on the distribution of pelagic species demonstrate that SST fronts are important habitat features (Sims and Quayle , Becker et al. , Scales et al. ).
Except for the large and established fronts, the North Sea fronts are very dynamic, and environmental conditions can change over relatively short temporal and fine spatial scales (Pingree et al. ). For our model it was important to capture the day‐to‐day movements of harbor porpoises foraging for patchy prey items, which in turn are responding to short‐term oceanographic variations. Longer term, seasonal effects could be captured effectively using SST and day length in the model by subsequently averaging the predicted daily harbor porpoise densities within each season. Physiographic or static predictors, such as water depth, slope or distance to coast, were good proxies for harbor porpoise distribution, also in combination with dynamic predictors like tidal range or currents as described for western Scotland and Wales (Pierpoint , Marubini et al. , Embling et al. , Isojunno et al. , Booth et al. ). In other areas tidally driven upwelling zones (Skov and Thomsen ), fronts (Johnston et al. ) and mean seasonal SST (Gilles et al. ) have also been shown to be linked to the distribution and movements of harbor porpoises.
Limitations and caveats
The SCANS II survey was conducted only during one summer, and no comparable survey data were available in the northern part of the study area during spring and fall. Our use of the multiseason data set to model porpoise–habitat relationships assumed that at least some of these relationships would be stable across seasons. This assumption cannot presently be evaluated, but we were appropriately cautious with our interpretation and have attempted to characterize the underlying uncertainty in our density maps. These spatially explicit estimates of uncertainty show that density predictions in regions and seasons with little or no effort have greater uncertainty (see Figs. 5–7). In addition, we explicitly accounted for complex sources of uncertainties throughout all steps of data analysis and model development.
The explained deviance of 24% is relatively good for cetacean habitat modeling studies (cf. e.g., Becker et al. , Mannocci et al. , Forney et al. ), although a considerable part of the variation in the porpoise distribution remained unexplained. It would be naïve to believe that we would be able to include all relevant drivers of porpoise distributions, or that we have necessarily identified the best proxy variables to use for model‐based predictions. Data on prey abundance may have increased explanatory power, but these data are currently not available at the appropriate spatial and temporal scales. Additionally, environmental dynamics that aggregate prey are often more consistent and Torres et al. () found in a direct comparison that their model was more successful when using environmental variables instead of prey data as predictor variables.
Other predictors we considered were either (1) not available for the complete survey period (e.g., wind stress (upwelling) from QuikSCAT, only available through November 2009, and its replacement ASCAT showed persistent differences during the short period of overlap), (2) had large data gaps on daily or even weekly resolution (e.g., chlorophyll‐a), or (3) were not appropriate for our study area (e.g., some studies include deviation of sea surface height (Palacios et al. , Forney et al. ); however, remotely sensed altimeter data are not useful for the North Sea due to the strong tides).
Management and conservation perspectives
This study currently presents the most comprehensive, seasonal habitat‐based estimates of porpoise density in the North Sea, and it is intended to support the management of anthropogenic impacts (e.g., regulate timing of pile‐driving for offshore renewable energy installations). This project demonstrates that annual small‐scale (i.e., EEZ‐wide) surveys are important to assess trends and seasonal distribution patterns of small cetaceans. However, as national waters do not represent biologically meaningful management units, a coordinated monitoring effort and the rigorous application of standardized methods is needed to fulfill requirements of the EU Marine Strategy Framework Directive. To allow assessment of broader‐scale population distribution shifts, a SCANS‐type survey with larger spatial coverage should be undertaken on a more regular basis (e.g., every 6–8 yr). This would also benefit the 6‐yr time scale that is implied in some EU legislation, for example, the Habitats Directive. Additionally, knowledge of species–habitat relationships is important for the assessment of stock status and trends, because it allows the explicit incorporation of habitat information when estimating abundance and identifying management units. In the context of marine spatial planning, a habitat‐based model is required since it will allow for predictions of porpoise densities within areas of special interest for any day within the original survey time span. Using day length as a proxy for seasonality makes such a model a valuable tool for any mesoscale‐planning endeavor within the study region.
A challenge for scientists and regulators is the assessment of cumulative effects of multiple anthropogenic stressors on population dynamics of marine mammals and the prediction of population‐level consequences of disturbance (e.g., Nowacek et al. ). This will become increasingly important in the face of a changing climate and increased anthropogenic development in marine habitats and hence to ensure a sustainable management of the population.
Future work
Forecasting is an emerging methodological application within the field of SDM and applies new combinations of consolidated methods. Advances in remote sensing and ocean model products allow for predictions of species density at a variety of temporal scales. This is important for the assessment of current and future risks to marine mammals in a changing environment. Becker et al. () showed that the Regional Ocean Modeling System (ROMS) SST forecasts were promising for forecasting of cetacean density. The advantage of our model is that it uses only SST as a dynamic predictor, thus, a similar approach could be feasible, for short‐term and by including climate change scenarios.
The season‐specific population density maps presented here are an essential input for individual‐based models (IBMs), where population dynamics can be predicted based on simulated animals that move in the landscape in a realistic manner (Grimm and Railsback ). For the harbor porpoise an IBM has already been developed for the inner Danish waters (Nabe‐Nielsen et al. , ), where local porpoise densities were used as a proxy for food availability. Animals were simulated to respond to disturbances by being displaced, resulting in reduced food intake and potentially reduced population sizes. The output of this study will make it possible to re‐parameterize this model for North Sea conditions. This will be an important step toward ensuring that harbor porpoises are managed sustainably in the presence of bycatch, noise from ships, wind‐farm constructions and other types of anthropogenic stressors, which notably overlap to a great degree in the North Sea (Halpern et al. ).
Acknowledgments
This study is part of the DEPONS project (
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
© 2016. This work is published under http://creativecommons.org/licenses/by/3.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Abstract
Effective species conservation and management requires information on species distribution patterns, which is challenging for highly mobile and cryptic species that may be subject to multiple anthropogenic stressors across international boundaries. Understanding species–habitat relationships can improve the assessment of trends and distribution by explicitly allowing high‐resolution data on habitats to inform abundance estimation and the identification of protected areas. In this study, we aggregated an unprecedented set of survey data of a marine top predator, the harbor porpoise (Phocoena phocoena), collected in the
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
Details
1 Institute for Terrestrial and Aquatic Wildlife Research, University of Veterinary Medicine Hannover Foundation, Büsum, Germany; Marine Mammal and Turtle Division, Southwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, La Jolla, California 92037, USA
2 Institute for Terrestrial and Aquatic Wildlife Research, University of Veterinary Medicine Hannover Foundation, Büsum, Germany
3 Marine Mammal and Turtle Division, Southwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Santa Cruz, California 95060, USA
4 IMARES Wageningen Institute for Marine Resource & Ecosystem Studies, Den Helder, The Netherlands
5 Royal Belgian Institute of Natural Sciences (RBINS), Operational Directorate Nature, Ostend, Belgium
6 Department of Bioscience, Aarhus University, Roskilde, Denmark
7 IMARES Wageningen Institute for Marine Resource & Ecosystem Studies, Den Helder, The Netherlands; Department of Aquatic Ecology and Water Quality Management, Wageningen, The Netherlands